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ABSTRACT 

A  three-component,  broad-band  seismic  system  which 
records  directly  onto  seven  track  digital  tape  has  been 
developed  in  the  Geophysics  electronic  laboratories  at  the 
University  of  Alberta.  A  differential  filtering  technique 
using  low  noise  operational  amplifiers  extends  the  usable 
velocity  sensitivity  of  a  Willmore  Mk  II  seismometer  having 
a  natural  period  of  1.5  seconds  out  to  about  40  seconds 
under  normal  operating  conditions.  A  wide  spectrum  of 
seismic  energy  may  thus  be  recorded. 

With  the  availability  of  three  component  records, 
digital  processing  techniques  which  utilize  the  polarization 
properties  of  elastic  body  and  surface  waves  may  be  applied 
to  the  data  for  enhancement  of  the  signal  to  noise  ratio. 

The  application  of  two  such  processors  is  considered.  The 
first,  particularly  suited  for  enhancement  of  shorter 
period  body  waves,  is  applied  in  the  time  domain;  while  the 
second  filter  is  used  to  separate  the  longer  period  surface 
Rayleigh  and  Love  waves  by  modification  of  the  discrete 
Fourier  coefficients  of  the  data  traces  in  the  frequency 
domain.  Results  obtained  from  the  seismograms  processed 
by  these  methods  show  that  a  very  good  separation  of  phases , 
on  the  basis  of  the  polarization  of  the  seismic  energy, 
can  be  accomplished.  Better  interpretation  of  both  source 
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characteristics  and  localized  recording  conditions  is 
thus  possible. 
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CHAPTER  1 

POLARIZATION  FILTERS  -  A  REVIEW 

1.1  Introduction 

As  any  seismograph  record  is  always  noise  con¬ 
taminated,  the  detection  and  interpretation  of  seismic  events 
requires  knowledge  of  the  characteristics  of  both  signal  and 
noise.  Elastic  body  waves,  which  may  be  generally  separated 
into  P  (compressional )  and  S  (shear)  phases,  can  be  con¬ 
sidered  as  non-dispersive  group  arrivals  with  maximum  power 
in  the  0.3  to  10  second  period  range.  Surface  Rayleigh  and 
Love  waves  may  be  described  as  dispersive  group  arrivals  with 
maximum  power  in  the  2  to  100  second  period  range  for  earth¬ 
quakes  of  moderate  magnitude,  the  observable  periods  extending 
up  to  57  minutes  for  larger  teleseismic  events. 

Superimposed  on  these  signals  is  microseismic  back¬ 
ground  noise  as  well  as  signal  generated  noise.  Signal 
generated  noise  is  the  result  of  multiple  reflections  and 
refractions  of  P  and  S  body  waves  at  crustal  interfaces  and 
inhomogeneities,  and  originates  primarily  under  the  recording 
station.  Microseismic  noise,  which  is  considered  to  consist 
mainly  of  fundamental  and  higher  mode  Rayleigh  waves,  has  been 
shown  to  exhibit  a  sharp  peak  in  the  5  to  8  second  period 
range  (Brune  and  Oliver  1959). 
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As  a  result  of  this  sharp  peak  in  the  microseismic 
noise  spectrum,  frequency  bandpass  filtering  is  often  employed 
to  improve  the  signal  to  noise  ratio.  Initial  processing  of 
seismograph  records  for  body  wave  studies  often  involves  a 
bandpass  filter  with  the  low  frequency  cutoff  at  3  to  5 
seconds  period,  while  a  similar  cutoff  is  applied  as  a  high 
frequency  corner  to  records  used  for  the  investigation  of 
longer  period  surface  waves.  If  the  seismograph  data  is 
available  in  digital  form,  zero  phase  shift  recursive  digital 
filters  are  most  often  used  for  this  handpass  filter 
processing . 


However,  although  filtering  of  this  type  is  very 
effective  in  removing  microseismic  background  from  both  long 
and  short  period  data,  it  often  cannot  distinguish  between 
signal  and  signal  generated  noise.  This  is  especially  true 
when  the  higher  frequency  body  waves  are  considered.  The 
reason  for  this  is  apparent  when  we  consider  that  signal 
generated  noise,  arising  from  reflected  and  refracted  P  and 
S  body  waves,  can  be  expected  to  exhibit  a  frequency  spectrum 
similar  to  the  original  phases. 

Although  bandpass  filtering  serves  to  improve  the 
quality  of  seismic  records  by  enhancing  the  signal  to  noise 
ratio,  difficulty  may  still  arise  in  the  attempted  identi¬ 
fication  of  phases  whose  frequency  characteristics  are  similar. 
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1.2  Polarization  of  Seismic  Energy 

A  further  property  of  seismic  events  which  may  be 
used  for  the  enhancement  of  recorded  data  is  that  of  signal 
polarization.  Both  compressional  and  shear  waves  exhibit  a 
high  degree  of  linear  polarization.  Particle  motion  coin¬ 
cides  with  the  direction  of  propagation  of  energy  for  com¬ 
pressional  (P)  phases  and  is  perpendicular  to  the  direction 
of  propagation  for  the  transverse  (S)  phases.  Surface  waves 
of  the  Rayleigh  type  are  generally  elliptically  polarized  in 
the  vertical  -  radial  plane,  the  fundamental  modes  dis¬ 
playing  retrograde  particle  motion  and  the  higher  modes  pro¬ 
grade  or  retrograde  ellipticity.  Surface  Love  waves  are 
also  found  to  be  rectilinear ly  polarized,  but  in  a  horizon¬ 
tal  plane  orthogonal  to  the  direction  of  wave  propagation. 

Since  most  microseismic  background  is  of  the 
Rayleigh  type,  it  will  exhibit  elliptical  polarization,  but 
with  little  preferred  directionality.  Signal  generated 
noise  may  also  be  polarized,  but  again  the  direction  of 
polarization  is  random  in  nature.  By  using  these  various 
characteristics  of  polarized  particle  trajectories,  filters 
may  be  designed  which  preserve  motion  if  it  satisfies 
specified  conditions  of  polarization  in  a  particular  direc¬ 
tion,  and  which  attenuate  motion  that  does  not  satisfy  the 


desired  criteria. 
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For  example,  a  filter  may  be  designed  so  that  it 
passes  only  energy  which  exhibits  a  high  degree  of  polariz¬ 
ation  (or  high  rectilinearity )  in  the  direction  of  incident 
compressional  pulses.  In  this  case,  P-type  motion  would  be 
preserved.  On  the  other  hand,  high  rectilinearity  in  a 
direction  perpendicular  to  the  direction  of  propagation  of 
the  seismic  waves  might  be  of  interest,  so  that  a  filter 
could  be  designed  to  pass  motion  which  is  SV  in  nature. 

Both  of  these  filters,  of  course,  would  attenuate  signal 
which  was  not  highly  rectilinear  in  the  specified  direction. 

We  see  then,  that  the  more  or  less  randomly  polar¬ 
ized  microseismic  background  and  signal  generated  noise  would 
be  attenuated,  as  well  as  the  elliptically  polarized  Rayleigh 
waves  and  the  transversely  polarized  Love  phases.  P  and  S 
body  wave  phases  would  be  preserved. 

For  applications  where  the  main  interest  might  be 
in  surface  Rayleigh  and  Love  type  arrivals,  similar  filters 
may  be  designed.  These  filters  do  not  require  the  high 
degree  of  linear  polarization  displayed  by  the  body  wave 
filters,  but  rather  weight  the  input  according  to  how  closely 
the  particle  motion  trajectory  corresponds  to  theoretical 
Love  and  Rayleigh  wave  motion  arriving  from  some  pre-assigned 
direction.  Simons  (1968)  for  example,  considers  rectiline¬ 
arity  of  Rayleigh  waves  based  on  a  theoretical  horizontal/ 
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vertical  displacement  ratio  of  about  0.8  (Dorman  and 
Prentiss ,  1960 ) . 

1.3  Polarization  Filters 

Considering  these  patterns  of  particle  motion, 
various  workers  have  designed  polarization  filters  for  the 
enhancement  of  seismic  data.  Shimshoni  and  Smith  (1964) 
suggest  that  the  time  averaged  cross  product  of  vertical  and 
radial  components  of  ground  motion  may  be  used  as  a  measure 
of  rectilinearity .  The  radial  trace  is  obtained  from  a 
rotation  of  the  two  horizontal  components  into  a  path  corres¬ 
ponding  to  the  great  circle  azimuth  from  source  to  receiver, 
and  the  computed  cross  product  is  multiplied  by  the  original 
signal.  A  function  of  ground  motion  which  enhances  recti- 
linearly  polarized  signal  is  thus  obtained. 

Another  process  described  by  the  same  workers 
decomposes  the  vertical  and  radial  motions  into  their  Fourier 
components  and  computes  the  parameters  of  an  equivalent  el¬ 
lipse  at  each  instant  in  time.  Eccentricity,  major  axis,  and 
angle  of  inclination  of  the  ellipse  from  the  vertical  are 
displayed  for  the  frequency  at  which  maximum  power  is  arri¬ 
ving.  This  information  is  used  to  provide  criterion  for  the 
identification  of  P  and  SV  type  motion.  In  both  processes 
described,  no  use  is  made  of  the  transverse  component  of 
ground  motion. 
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A  particular  polarization  filter  known  as  REMODE 
(REctilinear  Motion  DEtection)  which  has  received  considerable 
attention  was  originally  developed  by  Mims  and  Sax  (1965)  at 
Teledyne  Industries,  with  subsequent  work  by  Griffin  (1966a,  b) . 
This  process  uses  the  highly  rectilinear  properties  of  P  and 
SV  motion  to  discriminate  against  the  elliptically  or  randomly 
polarized  nature  of  signal  generated  noise.  Inputs  Z  (vertical) 
and  R  (radial)  are  rotated  so  that  the  expected  direction  of 
incident  body  waves  bisects  the  angle  between  the  two  ortho¬ 
gonal  components  (Figure  1.1).  In  this  way,  all  P  and  SV 
motion  is  represented  in  total  by  these  two  rotated  components 
so  that  a  similar  shape  of  the  Z  and  R  waveforms  (signal) , 
contrasts  with  dissimilar  shape  if  noise  is  present. 

An  estimate  of  these  similarities  is  obtained  by 
computing  the  cross-correlation  function  C(T),  of  Z (t)  and 
R(t),  over  a  specified  window  centered  at  some  time  t  on  the 
record.  If  the  polarization  in  the  R-Z  plane  is  predominan¬ 
tly  rectilinear,  the  even  part  of  C(T)  will  be  large  relative 
to  the  odd  part.  If  motion  of  elliptical  and  random  polar¬ 
ization  prevails,  the  even  part  of  CfT)  will  be  small  relative 
to  the  odd  part.  By  convolving  the  even  part  of  C(T)  with  the 
original  time  series  R(t)  and  Z (t) ,  motion  of  high  rectiline- 
arity  is  enhanced  and  elliptically  polarized  motion  is  atten¬ 


uated  . 
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Figure  1.1.  Rotation  of  Z  (vertical)  and  R  (radial)  components 


for  REMODE  processing  (after  Fuchs  1969) . 
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REMODE  as  described  above  will  pass  both  the  line¬ 
arly  polarized  P  and  SV  type  motion.  Fuchs  (1969)  and 
Basham  (1967)  describe  a  modification  originally  suggested 
by  Griffin,  which  serves  to  "tune"  the  filter  for  P  or  SV 
phases.  With  the  configuration  of  the  rotated  Z  and  R 
components  as  in  figure  1.1  we  find  that  for  P  or  compres- 
sional  motion, 

R (t) Z ( t) >0  (1.1) 

and  for  SV  motion  perpendicular  to  the  direction  of  energy 
propagation , 

R (t) Z (t) <0  (1.2) 

These  relationships  are  illustrated  in  figure  1.2. 


Thus,  if  compressional  motion  is  of  interest,  an 
operator  0  may  be  defined  such  that 

ir 


0p  =  1,  R ( t ) Z ( t) >0 
0p  =  0,  R ( t) Z (t) <0 


(1.3) 


and  applied  to  the  REMODE  filtered  traces  for  the  preser¬ 
vation  of  P-motion.  Conversely,  an  operator  0gv  defined  so 
that 


0Sv  =  1,  R ( t ) Z ( t) <0 
0SV  =  0,  R ( t ) Z (t) >0 


(1.4) 


may  be  applied  to  the  filtered  components  if  SV-type  motion 
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SIGNS  OF  COMPONENTS 
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Figure  1.2  Detection  of  P  and  SV  motion  by  the  product  of 
the  Z  and  R  components  of  motion  (after  Basham  1967). 
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is  desired. 


Archambeau  and  Flinn  (1965)  point  out  that  a  REMODE 
filter  of  the  type  described  above  can  be  approximated  in  the 
frequency  domain  by  a  cosine  function  of  the  form 

<Moo)-  cos[<J>R(aj)  -  <J>Z  (oj)  ]  (1.5) 

where  4>r(oj)  and  <j>z  (to)  are  the  phase  angles  for  the  frequency 
components  R(co)  and  Z  (co)  of  the  original  traces.  This  filter 
function  will  admit  displacements  in  phase  or  180°  out  of  phase 
on  the  vertical  and  radial  components  so  that  purely  recti¬ 
linear  motion  has  unit  gain  and  other  motion  is  attenuated. 
Elliptically  polarized  motion  however  may  still  be  passed,  as 
for  example  cf)R  (oo) -<J>Z  (oo)  =  60°  ,  in  which  case  0  (co)  =  0.5. 

Higher  order  REMODE  filters  with  responses  approximated  by 

$(co)-  cos11  [<J)R  (  go  ) -4>z  (oj)  ]  (1.6) 

have  also  been  investigated  (Fuchs  1969). 

Flinn  (1965)  describes  a  processor  for  three  com¬ 
ponent  seismic  records  which  considers  not  only  rectiline- 
arity,  as  does  REMODE  and  the  earlier  filter  described  by 
Shimshoni  and  Smith,  but  also  rectilinearity  in  some  parti¬ 
cular  direction  in  three  dimensional  space.  The  two  hori¬ 
zontal  components  of  ground  motion  are  rotated  so  that  the 
radial  component  corresponds  with  the  expected  great  circle 
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azimuth  from  source  to  receiver  and  the  transverse  component 
is  orthogonal  to  this.  Over  a  specified  time  window  of 
length  NAt,  where  At  is  the  digitizing  interval  of  the  data, 
parameters  of  an  equivalent  ellipsoid  are  calculated  in  the 
time  domain  from  the  covariance  matrix  for  the  set  of  N 
points  (Z-j_,  Rj_,  Tj_)  i  =  1***N.  Z  (vertical)  ,  R(radial)  and 
T (transverse)  represent  the  three  orthogonal  components  of 
ground  motion. 

The  ratio  of  the  lengths  of  the  major  and  secondary 
axes  of  the  ellipsoid  is  used  as  a  measure  of  rectilinearity 
and  the  orientation  in  space  of  the  major  axis  is  used  to 
determine  the  direction  of  polarization.  Appropriate  func¬ 
tions  of  these  two  properties  (rectilinearity  and  direction) 
are  then  applied  to  each  of  the  three  orthogonal  records. 

If  for  example,  rectilinearity  is  high  and  the 
major  axis  of  the  ellipsoid  is  found  to  lie  in  the  Z 
direction,  motion  on  this  component  will  be  preserved  when 
the  filtering  is  applied,  while  the  radial  and  transverse 
motions  will  be  attenuated.  Noise,  of  course,  since  it  can 
be  expected  to  display  neither  a  high  degree  of  polarization, 
nor  any  preferred  direction  will  be  attenuated  on  all  three 
records . 


Unlike  REMODE  or  the  process  described  by  Shimshoni 
and  Smith,  this  particular  polarization  filter  also  considers 
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the  transverse  component  of  ground  motion.  Although  the 
information  contained  by  the  transverse  record  is  not  as 
useful  as  that  contained  by  the  vertical  and  radial  compon¬ 
ents,  it  might  be  valuable  for  gaining  information  about 
local  conditions  near  the  recording  station.  This  could 
especially  be  the  case  if  local  conditions  give  rise  to 
large  horizontally  polarized  shear  phases  (SH)  due  to 
reflection  and  refraction  of  P  and  SV  body  waves. 

Polarization  filters  of  the  type  described  above 
can  be  more  useful  than  frequency  filters  for  enhancement 
of  signal  to  noise  ratio  particularly  when  the  signal  and 
noise  have  similar  spectral  characteristics.  It  must  be 
remembered  however,  that  these  types  of  filters  are  time- 
varying  nonlinear  processors  and  their  applications 
require  more  computer  time  than  for  example  a  digital 
frequency  filter  whose  response  does  not  change  with  time, 
and  which  may  be  applied  recursively  in  the  time  domain. 

1.4  Some  Applications  of  Polarization  Filters 

Various  workers  have  applied  polarization 
filtering  techniques  to  recorded  seismic  data  for  improve¬ 
ment  of  the  signal  to  noise  ratio.  Lewis  and  Meyer  (1968) 
apply  a  phase  filter  of  the  REMODE  type  as  described  by 
Archambeau  and  Flinn  (1965)  to  data  recorded  during  the 
Early  Rise  experiment  in  the  summer  of  1966.  The  filtering 
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is  applied  to  the  vertical  and  radial  components  of  ground 
motion  and  used  to  enhance  body  phases  at  both  the  compres- 
sional  and  shear  wave  portions  of  the  seismogram. 

Fourier  transforms  Fz  (oo)  and  FR(to)  of  the  vertical 
and  horizontal  motion  respectively  are  computed  and  used  to 
form  a  function  FZR(m)  =  Fz  (to)  FR(to)*,  (*denoting  the 
complex  conjugate  of  FR(to)).  The  argument  of  FZR(co)  will  be 
0  or  it  if  the  spectral  components  of  the  vertical  and  radial 
motion  are  in  phase  or  out  of  phase  and  will  be  tt/2  for 
Rayleigh  type  motion.  The  enhancement  of  body  phases  then  is 
obtained  by  considering  a  filter  function  Ff(to)  equal  to 
the  cosine  of  the  argument  of  FZR(oo).  That  is: 

Re[FZR(0J)  ] 

Ff(oo)  =  -  (1.7) 

ABS [Fzr(w) ] 


or 


Ff(oo)  =  coscf) 

<t>  =  ARG[FZr(oo)  ] 


(1.8) 


This  function  is  now  multiplied  by  the  original 
Fourier  transforms  so  that 


Fz  '  (to)  =  F z  (to)  -Ff  (to) 


Fr'  (to)  =  Fr  (to)  -Ff  (to) 


(1.9) 


.  ' 
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and  the  inverse  Fourier  transform  is  taken  to  obtain  the 
filtered  record  in  the  time  domain. 

In  the  application  by  Lewis  and  Meyer,  the  filter 

function  was  raised  to  the  seventh  power  and  applied  to  data 

digitized  at  100  samples/second.  A  1  or  2  second  window 

moved  along  the  record  in  steps  of  several  tenths  of  a  second 

* 

was  used  for  this  analysis. 

In  another  application,  Basham  and  Ellis  (1969) 
use  a  REMODE  filter  designated  as  a  P-Detection,  (P-D)  filter 
to  process  P-wave  codes  of  numerous  seismic  events  recorded 
in  western  Alberta.  Data  is  digitized  at  16.6  samples/second 
and  first  bandpass  filtered  0.25-2HZ  to  remove  microseismic 
and  cultural  noise.  The  polarization  filter  applied  to  25 
seconds  of  record  following  the  P  onset  allows  identification 
of  numerous  compressional  wave  arrivals  including  P,  pP  and 
sP .  P  P  and  PKP  phases  for  events  at  appropriate  epicentral 
distances  are  also  detected.  Determination  of  P-pP  times 
permits  focal  depth  calculations. 

Simons  (1968)  describes  a  process  which  examines 
three  component  particle  motion  trajectories  of  longer  period 
Love  and  Rayleigh  waves.  Apparent  horizontal  azimuth, 
eccentricity  and  phase  difference  between  the  vertical  and 
radial  components  is  computed  in  the  frequency  domain  from  the 
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discrete  Fourier  coefficients  of  the  data  traces  considered 
over  a  specified  time  window.  The  amplitude  coefficients 
are  then  modified  by  appropriate  functions  of  eccenticity, 
apparent  azimuth  and  vertical-radial  phase  difference 
according  to  how  closely  each  particular  component  satisfies 
expected  theoretical  patterns  for  Love  and  Rayleigh  wave 
motion.  Using  these  modified  amplitude  coefficients,  the 
time  series  is  reconstructed  so  that  the  filtered  trace  is 
obtained  over  the  original  time  window.  The  window  is  then 
advanced  along  the  record  and  the  calculations  are  repeated. 

Chapter  4  describes  this  process  in  more  detail. 
The  polarization  filter  described  by  Flinn  (1965)  is 
considered  in  Chapter  3.  Examples  of  records  processed  by 
these  methods  are  presented.  A  new  seismic  system  develo¬ 
ped  at  the  University  of  Alberta  which  records  directly  onto 
digital  tape,  is  described  in  Chapter  2.  Examples  of 
three  component  recordings  from  this  system  are  also  given. 
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CHAPTER  2 

A  WIDE-BAND  SEISMOGRAPH  SYSTEM 

2.1  Characteristics  of  Standard  Seismographs 

The  wide-band  spectral  characteristics  of  seismic 
arrivals  from  earthquakes  as  well  as  the  dominant  peak  in 
microseismic  background  level  for  the  5-8  second  period 
range  was  described  in  section  1.1.  A  seismograph  system 
which  has  broad-band  frequency  characteristics  and  which  at 
the  same  time  attenuates  energy  in  the  region  of  major 
microseisms  is  thus  most  desirable  for  any  seismological 
observatory . 

To  meet  these  criteria,  many  seismograph  stations 
employ  both  a  short  and  a  long  period  recording  system.  The 
instruments  are  designed  so  that  the  long  period  end  of  the 
short  period  system,  and  the  short  period  cut-off  of  the 
long  period  system  attenuate  the  major  microseismic  band. 

In  most  cases,  seismometer  output  is  recorded  on  photo¬ 
graphic  paper  by  means  of  a  light  spot  deflection  of  a 
galvanometer.  The  frequency  characteristics  of  such  systems 
are  dependent  entirely  on  the  natural  resonant  period  and 
damping  constants  of  both  the  seismometer  and  the  galvano¬ 
meter  . 


Although  the  short  period  instrumentation  is 


' 
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generally  found  to  be  adequate,  realization  of  mechanical 
constants  required  to  obtain  long  period  operation  of  galva¬ 
nometers  and  seismometers  make  it  difficult  to  extend  the 
long  period  sensitivity  much  higher  than  about  150  seconds. 
Such  low  frequency  characteristics  necessitate  the  use  of 
large  and  heavy  galvanometer  coils  requiring  great  amounts 
of  power  to  produce  a  given  deflection,  or  seismometers  with 
very  large  masses. 

2.2  Use  of  Feedback 

Electrical  feedback  may  be  used  to  alter  the 
characteristics  of  a  seismograph.  Many  methods  utilize  an 
electromagnetic  transducer  which  couples  the  seismometer 
frame  and  suspended  mass  (Tucker,  1958;  Willmore,  1961;  de 
Bremaecker  et  al,  1962;  Sutton  and  Lathan,  1964).  The  feed¬ 
back  current,  i,  flowing  in  the  transducer  coil  causes  a  mass 
motion  equivalent  to  a  ground  acceleration  gi/m  where  g  is 
the  transducer  constant  and  m  the  suspended  mass.  The  effect 
of  this  simulated  ground  acceleration  is  an  alteration  of  the 
characteristics  of  the  system,  the  extent  of  the  alteration 
being  dependent  upon  the  transfer  function  of  the  feedback 
loop  as  well  as  the  transfer  function  of  the  system  without 
feedback . 


Willmore  (1961)  describes  a  system  in  which  the 


natural  periods  of  the  seismometer  and  galvanometer  are  close 
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together  near  the  center  of  the  total  response  curve.  The 
desired  bandwidth  is  obtained  by  means  of  negative  electrical 
feedback  to  the  seismometer.  The  amount  of  this  feedback  is 
proportional  to  the  velocity  of  the  galvanometer  spot  deflec¬ 
tion.  A  seismograph  constructed  with  this  type  of  feedback 
using  a  seismometer  and  galvanometer  each  with  a  natural 
period  of  20  seconds  is  shown  to  have  an  almost  flat  velocity 
response  from  0.5  to  800  seconds. 

Russell  et  al  (1968)  eliminate  the  need  of  a  trans¬ 
ducer  as  input  terminals  for  the  application  of  negative  feed¬ 
back  by  incorporating  the  seismometer  in  a  balanced  Maxwell 
impedance  bridge.  Voltage  applied  to  the  bridge  input  termi¬ 
nals  simulates  a  ground  acceleration  and  allows  the  introduc¬ 
tion  of  a  feedback  loop.  This  applied  voltage  is  the  output 
of  the  seismometer  modified  by  the  transfer  function  of  the 
feedback  loop.  Such  feedback  may  be  used  to  control  individ¬ 
ually  the  mass,  spring  constant  and  damping  constants  of  the 
seismometer.  The  feedback  arrangement  is  shown  to  increase 
the  resonant  period  of  a  Willmore  Mark  I  seismometer  from  1 
to  nearly  5  seconds.  It  should  be  pointed  out  however,  that 
although  these  methods  may  be  effectively  used  to  alter  the 
transfer  characteristics  of  a  seismograph  system,  no  improve¬ 
ment  in  signal  to  noise  ratio  is  possible  through  the  use  of 


such  feedback. 
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2.3  Wideband  Response  By  Filtering 

Wideband  frequency  characteristics  may  also  be 
realized  by  filtering  the  seismometer  output.  Pomeroy  et  al 
(1969)  describe  a  system  in  which  the  motion  of  a  pendulum 
with  a  natural  period  of  34  seconds,  critically  damped,  is 
monitored  by  an  electromagnetic  (velocity)  transducer  which 
drives  a  galvanometer-phototube  amplifier  system.  The 
galvanometer  (also  critically  damped)  has  a  natural  period  of 
100  seconds.  The  amplified  output  is  filtered  and  recorded 
photographically  as  well  as  by  a  strip  chart  recorder.  Two 
different  filters  may  be  used. 

Comparison  is  made  with  records  of  an  event  recor¬ 
ded  by  a  conventional  electromagnetic  seismograph  at  the 
same  location.  The  mechanical  constants  of  this  system  are 
the  same,  but  no  amplification  is  used.  The  higher  gain 
system  is  shown  to  have  recorded  several  long  period  body 
phases  not  detected  by  the  conventional  seismograph.  Rigid 
environmental  control  however,  is  required  for  the  high  gain, 
long  period  operation  and  the  system  is  enclosed  in  an  air¬ 
tight  chamber  in  a  deep  mine . 

At  this  university,  an  amplifier  has  been  designed 
by  M.  D.  Burke  in  the  Geophysics  electronic  laboratories  for 
use  with  a  Willmore  Mark  II  seismometer.  The  usable  velocity 
sensitivity  of  the  seismometer  having  a  natural  period  of  1.5 
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seconds  and  damping  0.6  critical  is  extended  out  to  the  40 
second  period  range  under  normal  environmental  conditions. 

The  desired  system  response  is  obtained  by  the  use 
of  second  order  active  R-C  network  filters  following  a  design 
criterion  described  by  Sallen  and  Key  (1955) .  Employing  low 
noise  operational  amplifiers,  a  differential  filtering 
technique  is  used  to  separate  broad-band  and  long  period 
data,  amplify  the  long  period  data  to  compensate  for  the  drop 
in  sensitivity  due  to  the  seismometer,  and  then  recombine  the 
signals.  The  two  signals  may  be  summed  in  phase  or  180°  out 
of  phase  so  that  microseismic  enhancement  or  attenuation  can 
be  obtained. 

Figure  2.1  shows  the  measured  frequency  and  phase 
response  of  the  amplifier  itself.  The  dashed  line  indicates 
the  case  of  "in  phase"  summation  while  the  solid  curve 
displays  the  attenuation  available  in  the  region  of  peak 
microseismic  background  when  the  two  signals  are  summed  in 
opposition.  With  appropriate  choice  of  circuit  components 
this  summation  point  may  be  placed  at  any  particular  period, 
thus  adapting  the  amplifier  to  the  predominant  background 
present  at  one  particular  recording  site. 

In  figure  2.2  we  see  theoretical  asymptotes  for 
the  velocity  sensitivity  of  the  complete  system  (seismo¬ 
meter  and  amplifier) .  Comparison  is  made  with  the  calibra- 
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Figure  2.1  Tripartite  system  amplifier  frequency  and 
phase  response.  The  dotted  line  indicates 
the  case  of  'in  phase'  summations,  the  solid 
line  the  case  of  'out  of  phase'  summation 
for  microseismic  attenuation. 


PHASE  SHIFT,(DEG)  MAGNITUDE  (db) 
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Theoretical  asymptotes  for  the  tripartite 
seismometer  -  amplifier  combination. 


Figure  2 . 2 


VELOCITY  SENSITIVITY 
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tion  curves  of  the  Edmonton  seismograph  station  operated  as 
part  of  the  Dominion  Observatory  network.  Instrumentation  in 
this  case  consists  of  the  standard  seismometer-galvanometer 
combination  as  described  in  Section  2.1.  Natural  periods  of 
a  Willmore  Mark  II  seismometer  and  galvanometer  are  1  second 
and  0.24  seconds  respectively  for  short  period  recording, 
while  for  long  period  operation  a  Press-Ewing  seismometer 
having  a  natural  period  of  29  seconds  is  combined  with  a  95 
second  period  galvanometer.  Extension  of  the  Willmore 
characteristics  for  long  period  operation  is  exemplified  by 
this  diagram. 

Intersection  of  the  short  and  long  period  response 
curves  for  the  Edmonton  station  indicate  the  microseismic 
attenuation  obtained  at  about  7  seconds  period  using  the  two 
systems.  With  the  new  system,  similar  attenuation  as  shown 
in  figure  2.1  is  available  by  filtering  the  output  of  a 
single  seismometer. 

Appendix  I  presents  a  detailed  description  of  the 
amplifier  transfer  function  as  well  as  a  derivation  of  the 
theoretical  velocity  sensitivity  curve  shown  in  figure  2.2. 
Theoretical  response  curves  for  the  transfer  function  have 
been  calculated  using  an  APL  program  and  are  also  given. 

2.4  Digital  Recording  System 


A  field  operable  digital  recording  system  designed 


for  three  data  channels  has  also  been  developed  in  the  Geo¬ 
physics  electronic  laboratories  at  this  university,  for  use 
in  conjunction  with  the  analog  instrumentation  just  described. 
Figure  2.3  shows  a  block  diagram  of  the  system. 

A  successive  approximation  analog  to  digital  con¬ 
vertor  is  used  to  produce  an  11  bit  binary  word,  including 
sign,  from  the  input  voltage.  When  conversion  is  complete, 
the  binary  word  is  transferred  to  a  storage  register  and  an 
end  of  conversion  (E.O.C.)  pulse  is  used  to  switch  to  the 
next  data  channel.  The  next  conversion  is  made  while  the 
previous  word  is  being  written  onto  digital  tape.  After 
A-D  conversion  of  the  third  channel  signal,  an  end  of 
frame  pulse  (E.O.F.)  generates  a  flag  bit  which  is  written, 
as  the  12th  bit,  onto  tape  with  the  binary  word.  This 
flag  bit  is  used  for  channel  identification  when  the  tape 
is  processed. 

The  recording  is  made  on  a  seven  track  incremental 
transport.  Two  seven  bit  characters  are  written  for  each 
conversion.  The  first  character  consists  of  six  bits  of 
the  binary  word  (most  significant  bits,  including  sign)  and 
a  parity  bit.  The  last  five  bits  of  the  binary  word, 
channel  identification  and  parity  bit  are  written  as  the 
second  character.  Odd  parity  is  used. 


A  crystal  oscillator  whose  output  is  divided  down 
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Figure  2.3  Block  diagram  of  the  digital  recording 


system. 
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to  400  Hz  controls  pulse  switching  and  was  used  to  drive  the 
incremental  tape  transport  at  50  steps/second.  With  the  two 
character  word,  25  conversions/second  are  made  so  that  a 
sampling  rate  of  8 . 33/second/channel  is  obtained.  An  inter¬ 
record  gap  is  written  every  6144  conversions.  Write-time 
for  this  gap  is  short  enough  so  that  only  one  data  point  per 
channel  is  lost. 

2.5  Recorded  Data 

Three  components  of  ground  motion  (vertical,  east- 
west  horizontal  and  north-south  horizontal)  for  a  number  of 
seismic  events  were  obtained  while  the  system  described 
above  was  operated  over  a  period  of  several  months  at  the 
seismological  observatory  near  Edmonton.  Recording  at 
556BPI  with  a  Digi-Data  tape  transport,  the  system  could 
operate  unattended  for  approximately  3  days.  Helicorder 
records  of  the  vertical  component  served  as  an  analog 
monitor  for  determining  which  of  the  digitized  data  was  to 
be  used.  Initial  processing  consisted  of  an  unpacking 
procedure  which  transferred  the  sampled  data  from  seven  track 
field  tape  to  nine  track  digital  tape  for  use  with  an  IBM 
360/67  computer  at  the  University  of  Alberta  Computing  Center. 
A  Calcomp  plotter  output  of  the  three  component  records  was 
then  produced  so  that  a  visual  display  of  the  digitized  data 


was  obtained. 
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Figure  2.4  illustrates  such  an  event  from  a  magni¬ 
tude  5.9  earthquake  which  occurred  near  Taulaud  Island  in 
the  Philippines  on  January  30,  1969.  Universal  time  is  in¬ 
dicated  on  the  record.  An  epicentral  distance  of  103°  was 
calculated  using  CGS  preliminary  co-ordinates.  Arrival 
times  of  body  phases  as  predicted  by  the  Jef f reys-Bullen 
tables  were  determined.  A  well  defined  P-coda  is  indicated 
on  the  Z  component,  followed  by  good  quality  PP  and  PPP 
arrivals.  The  CGS  depth  determination  of  70  kilometers 
was  used  for  pP  time. 

Events  indicated  in  the  figure  as  1?  L  ,  P2  and  P3 
show  good  correlation  with  their  respective  pP ,  PP  and  PPP 
arrivals.  At  an  epicentral  distance  of  103°,  the  original 
compressional  pulse  undergoes  diffraction  at  the  earth's 
core.  The  PP  phase  since  it  does  not  penetrate  as  deeply, 
but  rather  is  reflected  once  at  the  surface,  is  not 
diffracted.  Scattering  of  energy  due  to  the  diffraction 
is  illustrated  by  the  initial  P  phase  and  its  coda  which 
exhibit  lower  amplitude  than  the  later  PP  phases  which 
traverse  the  crust  and  mantle. 

The  SKS-SKKS-S  group  of  waves  are  easily  identi¬ 
fied  on  the  North  and  East  horizontal  components. 

Multiply  reflected  PKKP ,  PKKS  and  a  particularly  clear 
SKKKS  arrival  on  the  East  component,  as  well  as  SS  and 
PPS  phases,  are  indicated  on  the  record.  The  onset  of 
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Figure  2.4  Digital  recording  of  a  magnitude  5.9  earthquake 
which  occurred  in  the  Philippines  on  January  30, 
1969.  A  =  103°.  Universal  time  is  indicated  in 
the  figure.  Various  events  have  been  identified 
on  the  basis  of  predicted  travel-times. 

At  =  0.12  seconds. 
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long  period  Love  and  Rayleigh  (LR)  waves  are  also  shown. 
Periods  range  from  about  35  -  40  seconds  at  the  beginning 
of  the  surface  wave  train  to  20  seconds  at  later  times. 

A  comparison  of  another  earthquake  detected  by 
this  system  with  the  same  event  recorded  by  the  standard 
seismometer-galvanometer  combination  onto  photographic 
paper  has  also  been  made.  In  this  case,  the  data  was  first 
processed  with  a  digital  zero  phase  shift  Butterworth 
filter  using  a  Fortran  IV  subroutine  developed  by  T. 

Alpasan  (1968)  at  this  university.  A  short  period  (0.25  - 
5  seconds)  and  a  long  period  (5  -  150  seconds)  pass-band 
have  been  used  in  order  to  match  as  closely  as  possible 
the  short  and  long  period  responses  of  the  conventional 
system.  Figures  2.5  and  2.6  show  the  impulse  and  amp¬ 
litude  response  of  these  two  filters.  In  actual  practice, 
the  gain  factor  introduced  by  the  filter  process  is 
removed  so  that  the  response  is  normalized  to  unity. 

Figure  2.7  illustrates  the  theoretical  system  response 
obtained  using  these  filters  and  compares  these  curves 
with  the  calibration  for  the  conventional  instrumentation. 

The  short  period  seismogram  from  the  standard 
system  and  the  0.25  -  5  seconds  filtered  record  from  the 
digital  system  is  shown  in  figure  2.8  (traces  marked  U,A.). 
The  event  is  a  magnitude  5.4  earthquake  at  an  epicentral 
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Figure  2.5  Impulse  and  amplitude  response  of  a  0. 2-4.0  Hz 
digital  zero  phase  shift  Butterworth  filter. 
The  impulse  response  has  been  computed  for 
512  points.  At=0.12  seconds. 
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Figure  2.6  Impulse  and  amplitude  response  of  a  0.007- 
0.2  Hz  (5-150  seconds)  digital  zero  phase 
shift  Butterworth  filter.  The  impulse 
response  has  been  computed  for  2048  points. 


At=0.12  seconds. 
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Figure  2.7  The  effect  of  the  digital  filters  shown  in  figures 
2.5  and  2.6  on  the  asymptotes  for  the  tripartite 
system.  The  solid  curve  marked  T.S.  is  for  the 
0.2-4  Hz  pass-band;  the  dashed  curve  for  the 
0.007-0.2  Hz  pass-band.  The  smooth  curves 
represent  the  calibration  of  the  short  period 
(Willmore)  and  long  period  (P-E)  systems. 
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distance  of  22°  which  occurred  in  the  Gulf  of  California  on 
March  21,  1969.  Predicted  arrival  times  of  P,  PP  and  PPP 
compressional  phases  for  a  surface-focus  earthquake  are 
indicated.  Since  the  depth  of  the  event  is  undetermined, 
these  times  may  be  2-5  seconds  late. 

The  5  -  150  seconds  filtered  record  is  compared 
with  the  long  period  seismogram  for  the  same  event  in 
figure  2.9.  To  aid  comparison,  the  original  vertical  (Z) 
component  of  the  conventional  system  is  shown  along  with 
tracings  of  the  three  records.  Predicted  arrival  times 
of  P,  S,  PCP  and  SCS  body  phases  are  indicated.  A  well 
defined  LR  (Rayleigh)  onset  is  shown  on  the  Z  and  North 
components  and  the  LQ  (Love)  surface  train  is  indicated 
on  the  East  recording. 

The  examples  presented  here  illustrate  the  wide¬ 
band  recording  capabilities  of  this  system.  In  figure  2.4 
we  see  on  the  same  record  a  well  defined,  short  period  P- 
coda  as  well  as  long  period  surface  Love  and  Rayleigh 
wave  trains.  Figures  2.8  and  2.9  illustrate  how  filtering 
may  be  applied  to  the  wide-band  digital  data  so  that  any 
desired  frequency  band  within  the  total  response  of  the 
system  may  be  extracted  from  the  original  record.  Al¬ 
though  this  system  does  not  extend  out  to  the  very  long 
periods,  the  wide-band  frequency  characteristics  and  good 
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Figure  2.8  Comparison  of  the  conventional  short  period 

seismogram  with  the  0.2-4  Hz  filtered  digital 
record  indicated  by  the  traces  marked  U.A. 
Three  components  of  ground  motion  are  shown. 
The  earthquake  of  magnitude  5.4  occurred 
in  the  Gulf  of  California  on  March  21,  1969 
at  04^56m20  SUT.  Universal  time  is  indicated. 


A=22° . 


GULFof  CALIFORNIA 
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Figure  2.9  Comparison  of  the  conventional  long  period 
seismogram  with  the  5-150  second  filtered 
digital  record  indicated  by  the  traces  marked 
U.A.  The  first  record  is  a  reproduction  of 
the  original  vertical  photographic  recording. 
The  second,  fourth  and  sixth  traces  are  copies 
of  the  three  component  photographic  seismograms 
presented  beside  the  U.A.  recordings  to  aid 
comparison.  The  event  is  the  same  Gulf  of 
California  earthquake  shown  in  figure  2.8. 
Universal  time  is  indicated. 


f 


GULF  of  CALIFORNIA  m  =  5.4,  march  21,  1969 
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dynamic  range  available,  along  with  its  adaptability  for 
operation  in  the  field  show  it  to  be  very  valuable  for  the 
detection  and  interpretation  of  seismic  events. 
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CHAPTER  3 

A  TIME  DOMAIN  POLARIZATION  FILTER 

3.1  Introduction 

Particle  motion  trajectories,  and  the  polarization 
characteristics  of  body  phases  and  surface  wave  groups 
have  been  discussed  in  Chapter  1.  A  number  of  time-varying 
nonlinear  processors  which  utilize  these  polarization 
properties  to  enhance  the  signal  to  noise  ratio  of  a  seis¬ 
mogram  have  been  described  (Shimshoni  and  Smith,  1964; 

Mims  and  Sax,  1965;  Flinn,  1965;  Fuch,  1969).  In  this 
chapter,  a  processor  which  can  be  particularly  useful  for 
enhancement  of  the  shorter  period  body  phases  is  discussed 
in  more  detail,  and  examples  of  records  processed  with  this 
filter  are  presented. 

Since  the  data  was  recorded  using  vertical,  east- 
west  and  north-south  horizontal  detectors  as  described  in 
the  previous  chapter,  the  two  horizontal  components  of 
ground  motion  were  first  rotated  so  that  a  radial  and  a 
transverse  record  were  obtained.  To  do  this,  the  expected 
great  circle  azimuth  from  source  to  receiver  was  computed 
and  a  transformation  applied  to  the  data.  If  cf>  is  the 
expected  azimuth  measured  clockwise  from  North,  then  the 
transformation  of  co-ordinates  is  given  by 
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R  =  E  sin  $  +  N  cos  cj) 

T  =  -E  cos  (J)  +  N  sin  <j> 


(3.1) 


where  E  and  N  represent  the  east-west  and  north-south  com¬ 
ponents  of  ground  motion  respectively.  Figure  3.1  illus¬ 
trates  this  transformation.  The  new  co-ordinate  system  is 
right-handed  so  that  R  x  T  =  Z,  with  Z  up  being  positive. 


3.2  The  Covariance  Matrix 


Adopting  the  method  described  by  Flinn  (1965),  a 
polarization  filter  that  considers  both  rectilinearity  and 
direction  of  particle  motion  has  been  developed.  In  order 
to  obtain  measures  of  these  two  quantities,  the  covariance 
matrix  for  a  set  of  N  points  taken  over  each  of  the  three 
orthogonal  components  of  ground  motion,  R,  T  and  Z,  is 
computed.  A  specified  time  window  of  length  NAt,  where 
At  is  the  sampling  interval,  is  thus  considered.  To 
determine  the  covariance  matrix  for  this  set  of  observa¬ 
tions,  the  means,  variances  and  covariances  must  be  calcu¬ 
lated  for  the  three  variables  R,  T  and  Z. 

We  define  the  mean  or  expected  value  of  N  obser¬ 
vations  of  the  random  variable  W  as 


N 

yw  =  1_  £  w-j_  =  E  (W) 

N  i=l 


(3.2) 
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Figure  3.1  Illustration  of  the  transformation  of  seismometer 
co-ordinates  so  that  the  east-west  and  north-south 
components  of  ground  motion  are  rotated  into  a  radial- 
transverse  configuration.  The  expected  great  circle 
azimuth  is  the  angle  cf>  measured  clockwise  from  North. 
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The  covariance  between  N  observations  of  two  variables  x1 
and  x2  is  given  by 

Cov [X ! ,  X2]  =  E [  (X1-E(X1))(X2-E(X2))]  (3.3) 

which  by  equation  3.2  can  be  written 

1  N 

Cov  [X  ,  X  ]  =  _  Z  (x.  -  y  x )  (x  .  -  y  )  (3.4) 

1  z  n  21  2 
N  1=1 

It  is  evident  that  Cov[Xj  ,  X2]  =  Cov[X2,  X 2 ]  .  The  quantity 
Cov[Xlf  Xx]  is  defined  as  the  variance  of  x2  and  by  equation 
3.4  may  be  written  as 

1  N  2 

Var[Xj]  =  -  Z  (xj-Mj)  (3.5) 

N  i=i 

The  matrix  with  Cov[X.y,  Xs]  in  its  row  and  s ^  column 

(Y,  s  =  1,  2,  •••  n)  is  the  covariance  matrix  for  the  set 
of  n  random  variables  X j ,  j  =  1 ,  2,  •••  n.  If  X  is  the 

vector  of  the  random  variables  and  y  the  vector  of  means 
for  each  of  these  variables,  the  covariance  matrix  V  is 
defined  by  (Jenkins  and  Watts,  Chapter  3) 

V  =  E[  (X-y) (X-y) t]  (3.6) 

The  superscript  t  indicates  the  column  transpose  of  the 


vector . 


This  can  be  written  as 


* 
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Var  [X1  ]  Cov[XlfX2]  •••  Cov[X  ,X  ] 
Cov[X2/X1]  Var[X2]  Cov[X2,Xn] 


(3.7) 


C°v[xn/xi]  Cov[Xn,X2]  •••  Var[xn] 


where  the  covariances  and  variances  are  defined  in  equations 
3.4  and  3.5.  As  Cov[X^,Xj]  =  Cov[Xj,X^]  this  matrix  is 
symmetric  about  the  main  diagonal.  For  our  case  of  three 
variables  R,  T  and  Z  considered  over  the  time  window  NAt, 
equation  3.7  is  written 


V  = 


Var [R]  Cov [  R ,  T  ]  Cov[R,Z] 
Cov [  R ,  T  ]  Var [T]  Cov[T,Z] 
Cov [R,  Z  ]  Cov [T , Z ]  Var [ Z ] 


(3.8) 


3.3  Measures  of  Rectilinearity  and  Direction 


If  the  covariance  matrix  given  by  equation  3.8  is 
diagonalized,  an  estimate  of  the  rectilinearity  of  particle 
motion  trajectory  over  the  specified  time  window  can  be 
obtained  from  the  ratio  of  the  principal  axes  of  this 
matrix.  The  direction  of  polarization  may  be  measured  by 
considering  the  eigenvector  of  the  largest  principal  axis. 
Since  the  eigenvalues  of  V  will  lie  along  the  principal 
axes ,  by  computing  the  eigenvalues  and  their  corresponding 
eigenvectors,  the  appropriate  functions  of  rectilinearity 
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and  directionality  can  be  determined.  To  this  end  it  was 
considered  that  if  A1  is  the  largest  eigenvalue,  and  A2 
the  next  largest  eigenvalue  of  the  covariance  matrix,  then 
a  function  of  the  form 

F(X,  ,A  )  =  I-,!!)"  (3.9) 

V 

would  be  close  to  unity  when  rectilinearity  is  high 

and  close  to  zero  when  the  two  principal  axes 
approach  one  another  in  magnitude  (low  rectilinearity) . 

The  direction  of  polarization  can  be  determined  by  con¬ 
sidering  the  components  of  the  eigenvector  associated  with 
the  largest  eigenvalue  with  respect  to  the  co-ordinate 
directions  R,  T  and  Z. 

To  illustrate  this,  figures  3.2a-d  show  some 
computations  applied  to  sets  of  data  in  two  dimensions. 
These  points  were  computed  for  an  ellipse  and  then  per¬ 
turbed  by  the  addition  of  random  noise.  In  3.2b  and  3. 2d, 
the  data  points  were  generated  from  a  45°  rotation  of  the 
original  ellipse  determined  for  3.2a  and  3.2b  respectively. 
The  computed  covariance  matrix,  correlation  coefficient 
p,  F  (X 1  ,  A  2 )  for  n  =  1,  and  the  eigenvector  of  the 
principal  axis,  E,  is  shown  for  each  case.  The  corre¬ 
lation  coefficient  is  determined  from  the  usual  definition 
(Jenkins  and  Watts,  Chapter  3) 
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2  (Cov [ X j  ,X  ]  ) 

Pl2  =  - - -  (3.10) 

Var [X} ] Var [X2 ] 

Comparing  3.2a  and  3.2b  we  see  that  in  the  first  case, 
Cov[X1,X2]  is  small  with  respect  to  the  variance  terms  on 
the  main  diagonal,  and  the  eigenvector  E  indicates  a 
preferred  direction  along  the  X2-axis.  In  the  second  case, 
Cov[Xx,X2]  is  larger  and  E  =  (0.726,  0.689)  shows  no 
preferred  direction  along  either  the  Xx  or  X2  co-ordinate 
axis.  In  both  instances,  F(X1,A2)  is  approximately  the 
same.  Similar  properties  of  the  covariance  matrix  and  the 
eigenvector  E  are  indicated  by  figures  3.2c  and  3. 2d,  but 
here  the  rectilinearity ,  F(A  ,X2),  is  low  in  both  cases. 

We  see  then  that  when  the  off-diagonal  terms  of  the 
covariance  matrix  are  significant,  the  diagonalization 
will  introduce  a  rotation.  After  rotation,  the  orienta¬ 
tion  in  space  of  the  principal  axis  of  the  covariance  matrix 
will  be  given  by  the  components  of  its  eigenvector  relative 
to  the  original  co-ordinate  system.  F(AlfA  )  will  give  an 
estimate  of  the  degree  of  polarization  along  the  major  axis. 
In  each  of  figures  3.2a-d,  the  data  points  could  be  con¬ 
sidered  as  representing  particle  motion  trajectory  over  a 
specified  time  in  one  of  the  orthogonal  planes  of  the  R,  T, 

Z  co-ordinate  system. 

Suppose  now  that  in  the  matrix  of  equation  3.8, 
the  covariance  terms  are  small  with  respect  to  those  on  the 


« 
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Figure  3.2  Computations  of  the  covariance  matrix,  correla¬ 
tion  coefficient  p,  F(XlfA  )  for  n=l,  and  the 
eigenvector  E  of  the  principal  axis  of  the 
covariance  matrix  for  sets  of  20  points  in  two 
dimensions.  The  data  was  generated  by  the 
addition  of  random  noise  to  a  perfect  ellipse. 
Figures  3.2b  and  3.2c  represent  rotations  of 
the  original  ellipse  by  45° . 


F  (X  x  ,X2)  =  0.950 
E  =  (-0.009,0.999) 


3.2b 


X 


2 


13.775  11.756 

[11.756  12.519 

p  =  0.895 

F(X1  ,\2)  =  0.945 

E  =  (0.726,0.689) 


2.416  0.221 

0.221  3.426 


p  =  0.077 
F(XirX2)  =  0.318 


2.803  0.420 

0.420  3.085 


p  =  0.143 
F  ( A 1  ,X2)  =  0.262 
E  = 


E 


(0.205,0.979) 


(0.584,0.812) 
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main  diagonal,  on  which,  Var [R] >Var [T] >Var [Z ] .  The  two 
principal  axes  would  thus  correspond  very  closely  to  Var[R] 
and  Var[T] ,  and  the  eigenvector  associated  with  the  major 
axis,  Var[R],  would  have  its  largest  component  in  the  R 
co-ordinate  direction.  If  Var [R] >>Var [T] ,  then  F(A1,X2) 
in  equation  3.9  would  be  close  to  unity  and  we  would  have 
high  rectilinearity  in  the  R  direction.  If  Var [R] War [T] , 
then  the  direction  of  polarization  would  still  be  pre¬ 
dominantly  along  the  R  co-ordinate  axis,  but  the  recti¬ 
linearity  would  be  low.  If  the  off-diagonal  terms  of  the 
matrix  were  significant,  the  diagonalization  would  introduce 
a  rotation,  and  the  orientation  in  space  of  the  major  axis 
would  be  given  by  the  components  of  its  eigenvector  relative 
to  the  original  R,  T,  Z  co-ordinate  system.  Depending  on 
the  amount  of  rotation  necessary  to  produce  this  diagonal¬ 
ization,  the  rectilinearity  would  exhibit  either  a  preferred 
direction  or  no  preferred  direction.  By  combining  F(A  ,X2) 
and  the  appropriate  eigenvector,  we  see  that  measures  of 
rectilinearity  and  directionality  can  be  obtained.  Deter¬ 
mination  of  a  suitable  filter  function  is  then  possible. 

3.4  The  Polarization  Filter 

Applying  this  analysis  to  the  digital  seismograms, 
the  covariance  matrix  is  computed  for  a  specified  time 
window  of  length  NAt  centered  about  tQ,  where  tQ  is  allowed 
to  range  over  the  entire  record  length  of  interest. 
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Eigenvalues  and  the  corresponding  eigenvectors  of  this 
matrix  are  then  determined.  The  measure  of  rectilinearity 
for  the  time  t  is  given  by 

RL(t0)  =  [F  ( A  t , A  2 ) ] J  (3.11) 

where  F(A  ,A2)  is  defined  in  equation  3.9.  If  we  represent 
the  eigenvector  of  the  principal  axis  with  respect  to  the 
R,  T ,  Z  co-ordinate  system  by  E  =  (elfe2,e3)  then  the 
direction  functions  at  time  tG  are  given  by 

Dr (t0)  =  (e,)K 

DT(t0)  =  (e2 ) L  (3.12) 

Dz(t0)  =  (e3)M 

Since  the  eigenvector  is  normalized  (|e|  =  1) ,  we  see  that 

for  each  of  these  functions  D. 

i 

0  <  Di  <  1  (3.13) 

The  exponents  J,  K,  L  and  M,  as  well  as  n  (equation  3.9) , 
are  determined  empirically.  Once  these  functions  are 
computed,  the  window  is  moved  down  the  record  one  sample 
interval  and  the  calculations  are  repeated.  When  the  end 
of  the  desired  record  is  reached,  these  quantities  are  then 
averaged  over  a  window  equal  to  about  half  the  original 
window  length.  This  has  the  effect  of  "smoothing"  these 
operators  so  that  contributions  due  to  any  anomolous  spikes 
are  subdued.  If  this  time  window  consists  of  M  points,  the 
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gain  functions  are  given  by 


L 


(3.14) 


M-l 


The  resulting  operators  are  then  used  as  a  point  by  point 
gain  control  to  modulate  the  rotated  records  so  that  at  any 
time  t,  the  filtered  seismograms  are  given  by 


Rf (t)  =  R ( t) -RL* (t) • Dr* (t) 


Tf (t)  =  T ( t )  • RL  *  ( t )  -DT*(t) 


(  3 . 15  ) 


Zf (t)  =  Z (t) ‘RL* (t) *DZ* (t) 


In  most  cases,  the  data  was  band-pass  filtered  with  a  zero 
phase  shift  digital  filter  before  rotation  so  that  the 
spectral  content  of  the  record  was  known.  The  window  length 
could  be  specified  so  as  to  be  consistent  with  one  or  two 
cycles  of  the  dominant  period.  Values  of  J,  K,  L  and  M 
which  appeared  quite  adequate  were  1,2,2 ,2  respectively. 

For  rectilinearity ,  n  =  0 . 5  or  1  was  used. 

3.5  Application  of  the  Polarization  Filter 


A  Fortran  IV  program  has  been  written  for  this 
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polarization  filter.  The  program  is  arranged  so  that  window 
length  for  both  the  covariance  matrix  (equation  3.8)  and  the 
smoothing  process  (equations  3.14),  starting  point  of 
calculations,  and  the  exponents  J,  K,  L,  M  and  n  can  be 
specified  as  input  data.  The  option  of  frequency  filtering 
is  also  available.  Figures  3.3  and  3.4  show  two  events 
processed  by  this  method.  The  first  example  is  of  a  mag¬ 
nitude  5.9  earthquake  which  occurred  in  the  Philippines  on 
January  30,  1969  at  a  depth  of  70  kilometers,  while  the 
second  is  the  nuclear  test  shot  'Benham'  of  magnitude  6.3, 
detonated  in  Nevada  on  December  19,  1968.  The  first  three 
traces  represent  the  vertical,  radial  and  horizontal 
components  of  ground  motion  after  rotation  (equation  3.1), 
followed  by  the  filter  functions  as  given  by  equations  3.14. 
The  last  three  represent  the  records  after  applying  these 
operators  as  in  equations  3.15.  The  filter  functions  and 
the  processed  records  are  displayed  from  the  point  in  time 
where  calculations  began. 

Data  in  figure  3.3  was  band-pass  filtered  0.3-4Hz 
before  rotation.  We  see  a  good  separation  of  events 
labeled  P l ,  P2  and  P3  in  the  P,PP  and  PPP  codas  of  the 
processed  seismogram.  The  good  correlation  of  these  events 
with  their  respective  pP  times  suggests  a  multiplicity  of 
sources  at  the  same  location  rather  than  a  single  release 
of  energy.  Attenuation  of  the  radial  and  transverse  traces. 
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Figure  3.3  Example  of  a  magnitude  5.9  earthquake  which 

occurred  in  the  Philippines  on  January  30,  1969, 
processed  by  the  polarization  filter  with  J=l, 
K,L,M=2  and  F  ( A  t  ,  A  2  )  =1-  (Ai.)  *5  .  The  event 

h  A  i 

occurred  at  10  29  40  UT .  Universal  time  is 
indicated  along  the  top  of  the  figure.  For  this 
event,  A=103° .  The  first  three  traces  are  the 
vertical,  radial  and  horizontal  components  of 
ground  motion  respectively;  the  next  four 
represent  the  computed  filter  operators;  the 
last  three  are  the  filtered  seismograms.  Data 
was  band-pass  filtered  0.3-4  Hz  before  rotation. 
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especially  during  the  P  and  PP  events  illustrates  how  the 
filter  enhances  motion  which  exhibits  a  preferred  direction 
of  polarization,  in  this  case,  the  Z  direction.  We  note 
also  the  increase  in  signal  on  the  two  horizontal  components 
and  a  corresponding  decrease  in  the  vertical  trace  at  the 
beginning  of  the  SKS-SKKS-S  group  of  waves  which  has  been 
identified  using  predicted  times  from  the  Jef f reys-Bullen 
travel-time  curves.  This  is  to  be  expected  since  little  or 
no  signal  should  be  present  on  the  vertical  component  of 
ground  motion  when  shear  (S)  waves  are  predominant. 

The  nuclear  test  shot  of  figure  3.4  was  recorded 
with  microseismic  attenuation  as  described  in  section  2.3. 
The  effect  of  this  is  seen  in  the  almost  zero  level  signal 
before  the  beginning  of  the  event.  In  this  case,  no  digital 
filtering  has  been  applied  to  the  data  before  rotation.  We 
see  here  that  the  processed  records  show  a  very  simple  P- 
coda,  as  would  be  expected  from  an  explosive  point  source. 
Again,  the  shear  phases  (S,SS,SSS)  are  present  only  on  the 
horizontal  components.  Of  interest  is  the  high  level  of 
activity  exhibited  by  the  radial  and  transverse  records 
near  the  beginning  of  the  event.  These  arrivals  cannot 
be  correlated  with  any  predicted  travel-times,  but  could 
possibly  be  due  to  body  phase  conversions  near  the  recording 
station . 


Figure  3.5  shows  a  particle  motion  diagram  for  the 
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Figure  3 . 


4  The  nuclear  test  shot  'Benham'  detonated  in 
Nevada  at  16h30mUT  on  December  19 ,  1968, 
processed  by  the  polarization  filter  with 
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J=l,  K,L,M=2  and  F(XlfX  )=!-( — ) 


Univer- 
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sal  time  is  indicated  along  the  top  of  the 
figure.  The  event  was  of  magnitude  6.3, 
with  A=16°.  The  display  of  traces  is  the 
same  as  in  figure  3.3.  No  band-pass  filtering 
has  been  applied  to  the  data. 
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vertical  (R-Z)  and  horizontal  (R-T)  planes  of  the  first  20 
seconds  of  this  event.  The  unfiltered  portion  of  the  three 
components  of  ground  motion  are  shown  with  time  in  seconds 
from  the  onset  indicated  above  each.  A  2.5  second  (21 
point)  window  is  represented  by  every  diagram.  This  window 
is  advanced  in  time  2.5  seconds  for  each  plot  so  that  a 
continuous  representation  of  the  event  is  displayed.  The 
build-up  of  energy  on  the  radial  and  transverse  records 
indicated  by  the  processed  seismograms  of  figure  3.4,  can 
be  seen  in  this  diagram  beginning  about  10  seconds  after 
the  onset  and  extending  out  to  about  20  seconds. 

The  R-Z  diagrams  for  the  onset  of  the  event  show 
that  the  original  compressional  pulse  is  arriving  at  near 
vertical  incidence,  but  the  orientation  of  the  principal 
axis  of  the  particle  trajectory  in  the  R-T  plane  suggests 
a  significant  departure  from  the  expected  great  circle 
azimuth  (R+  direction) .  Increase  in  the  transverse  com¬ 
ponent  from  10-15  seconds  indicates  that  the  horizontal 
motion  is  not  purely  SV  in  nature,  but  contains  a  signifi¬ 
cant  amount  of  horizontally  polarized  SH  motion  which  is 
being  detected  by  both  the  radial  and  transverse  seis¬ 
mometers.  Similar  effects,  particularly  at  this  recording 
location,  have  been  observed  by  Ellis  and  Basham  (1968)  and 
Basham  and  Ellis  (1969)  in  a  study  of  seismograms  taken  at 
various  stations  in  central  Alberta,  and  have  been 
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Figure  3.5  Particle  motion  diagrams  for  the  first  22.5 
seconds  of  the  nuclear  test  shot  'Benham'. 

The  unfiltered  portions  of  the  vertical, 
radial  and  transverse  records  are  shown  with 
time  in  seconds  from  the  onset  of  the  event 
indicated  above  each.  Sign  convention  is 
indicated  for  the  R-Z  and  R-T  planes  in  the 
first  diagram.  A  21  point  (2.5  seconds) 
window  is  represented  by  each  diagram  with  "1" 
representing  the  first  point  and  "21"  the  last 
point  so  that  the  direction  of  motion  can  be 
seen.  The  beginning  of  each  diagram  is 
separated  by  one  sample  interval  from  the  last 
point  of  the  preceeding  plot.  Time  for  point 
"1"  is  indicated  along  the  diagrams.  As  in 
figure  3.4,  no  band-pass  filtering  has  been 
applied  to  the  data  before  rotation. 
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suggested  as  being  the  result  of  body  wave  conversions  at 
or  near  the  base  of  the  sediments. 

In  the  course  of  the  present  study,  it  was  apparent 
that  a  larger  increase  in  transverse  energy  following  the 
original  compressional  pulse  was  associated  with  events 
arriving  at  azimuths  between  about  180-270  degrees.  This 
suggests  that  the  observed  transverse  motion  might  be 
generated  by  body  phase  conversions  at  major  topological 
features  in  the  vicinity  of  the  recording  station,  rather 
than  at  the  base  of  the  sediments.  Analysis  of  events 
arriving  from  various  azimuths  using  an  array  arranged 
around  the  present  recording  site  would  permit  a  more  reli¬ 
able  interpretation  of  these  observations. 

Although  the  particle  motion  diagram  of  figure  3.5 
shows  that  the  original  compressional  pulse  is  arriving  at 
near  vertical  incidence,  this  of  course  may  not  necessarily 
be  the  case  for  all  events.  If  the  angle  of  incidence  were 
inclined  appreciably  from  the  vertical,  then  pure  P  or  SV 
motion  would  be  present  on  both  the  vertical  and  radial 
components ,  and  separation  of  the  two  types  would  not  be  as 
efficient.  By  considering  a  unit  vector  A  =  (a1,a2,a3), 
which  is  orientated  at  the  observed  inclination  from  the 
vertical  in  the  R,T,Z  co-ordinate  system,  direction  functions 


of  the  form 
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where  E  =  (elfe2,e  )  is  the  eigenvector  of  the  principal 
axis  of  the  covariance  matrix,  could  be  obtained.  Rotation 
of  the  Z  and  R  components  of  motion  about  the  transverse 
horizontal  direction  so  that  the  former  vertical  instrument 
lies  in  the  P  direction  and  the  former  radial  instrument  in 
the  SV  direction  could  be  accomplished.  Application  of  the 
direction  operators  3.16  might  then  more  effectively 
separate  pure  P  and  SV  motion.  Transverse  motion  on  the 
original  radial  trace  would  also  be  eliminated  by  virtue  of 
the  rotation  in  the  vertical-radial  plane.  One  method  of 
determining  the  apparent  direction  of  incidence  is  by 
fitting  a  straight  line  to  the  vertical-radial  particle 
motion  considered  over  the  first  few  seconds  of  the  P-coda 
(Fuchs  1969 ) . 

A  comparison  of  filtered  and  unfiltered  records  from 
three  events,  A,B,C,  all  within  11  kilometers  of  each  other/ 
is  shown  in  figures  3.6  and  3.7.  The  earthquakes,  of 
magnitude  5.3,  5.4  and  5.5  respectively  occurred  in  the  Gulf 
of  California  within  about  a  3  hour  period  on  March  21,  1969. 
In  each  case,  the  data  was  band-pass  filtered  0.3-4HZ  before 
rotation.  On  the  seismograms  processed  with  the  polarization 
filter,  a  good  correlation  from  record  to  record  is  evident 
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Figure  3.6  Vertical,  radial  and  transverse  components  of 

ground  motion  for  three  events  which  occurred  in 
the  Gulf  of  California  on  March  21,  1969. 
Magnitudes  of  each  event  A,  B  and  C  are  5.3, 

5.4  and  5.5  respectively.  For  each,  A=22°. 
Universal  time  is  indicated  on  each  record. 
Digital  band-pass  filtering  of  0.3-4  Hz  has 
been  applied  before  rotation. 
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Figure  3 . 


7  The  seismograms  of  figure  3.6 
with  the  polarization  filter, 
events  from  predicted  P  and  S 
indicated.  For  these  records 
and  F(X1  ,X2)=  1- (^-) . 


after  processing 
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travel  times  is 
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for  the  P-PP-PPP  and  S-SS-SSS  wave  groups.  Since  these 
events  have  been  identified  from  predicted  times  for  a 
surf ace-f ocus  earthquake,  and  the  depth  of  each  is  unde¬ 
termined,  at  least  part  of  the  variations  in  arrival  time 
could  be  due  to  different  focal  depths.  The  increase  in 
horizontally  polarized  shear  motion  is  clearly  displayed 
at  the  beginning  of  the  Love  (LQ)  wave  train  on  the 
transverse  component  of  each  filtered  seismogram.  The 
onset  of  Rayleigh  (LR)  waves  is  also  clearly  indicated  on 
each  of  the  radial  components ,  and  a  separation  of  domi¬ 
nant  Love  and  Rayleigh  type  motion  can  be  seen.  The 
absence  of  any  vertical  (Z)  motion  over  the  portion  of  the 
filtered  seismogram  following  the  arrival  of  LR  indicates 
that  the  major  axis  of  polarization  never  lies  predomi¬ 
nantly  in  this  direction.  Assuming  the  Gutenberg  earth 
model,  the  theoretical  horizontal/vertical  displacement 
ratio  for  fundamental  long  period  Rayleigh  waves  is  0.8 
(Dorman  and  Prentiss,  1960).  We  would  thus  expect  the 
vertical  ground  motion  to  be  dominant  for  such  arrivals. 
Since  this  does  not  appear  to  be  the  case,  indications  are 
that  the  major  axis  of  the  particle  trajectory  for  the 
fundamental  Rayleigh  waves  does  not  lie  in  the  Z  direction, 
but  rather  is  inclined  at  some  angle  to  the  vertical. 
Haskell  (1953)  shows  that  transmission  of  Rayleigh  waves 
in  imperfect  elastic  medium  can  give  rise  to  such  effects. 
Whether  or  not  this  is  a  localized  condition  could  only  be 
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determined  by  comparison  with  seismograms  from  other 
recording  sites. 

The  examples  shown  here  illustrate  how  this  type 
of  polarization  filter  can  be  applied  to  enhance  motion 
which  is  rectilinear  in  a  preferred  direction.  Future 
applications  of  this  processor  to  data  sets  from  an  array  of 
recording  instruments  could  be  valuable  in  interpretation  of 
both  source  characteristics  and  localized  conditions  near 
the  receiver.  It  is  evident  however,  that  particle  trajec¬ 
tories  of  arrivals  such  as  Rayleigh  surface  waves  which 
exhibit  motion  more  elliptical  in  nature,  will  be  subject  to 
attenuation  along  one  of  the  major  axes.  This  is  seen  to 
be  the  case  in  figure  3.7  when  strong  Love  waves  exhibiting 
high  rectilinearity  in  the  transverse  direction  arrive  at 
about  the  same  time  as  the  Rayleigh  phases.  The  filter 
thus  appears  to  be  most  useful  for  the  enhancement  of  P  and 
S  body  waves . 
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CHAPTER  4 

A  FREQUENCY  DOMAIN  SURFACE  WAVE  DISCRIMINATION 

FILTER 

4 . 1  Introduction 

The  arrival  of  Love  and  Rayleigh  surface  waves  from 
a  teleseismic  event  is  generally  characterized  on  a  seismo¬ 
gram  by  an  increase  in  long  period  signal.  Since  the  group 
velocity  of  the  fundamental  Rayleigh  mode  is  less  than  that 
of  the  Love  mode  (Grant  and  West,  Chapter  3),  the  Rayleigh 
waves  will  generally  arrive  behind  the  Love  waves.  Most 
often  however,  this  time  difference  is  not  sufficient  to 
separate  the  Rayleigh  and  Love  groups  so  that  at  the  same 
time,  motion  which  is  elliptical  in  nature  will  be  present 
in  the  vertical-radial  plane  due  to  the  Rayleigh  waves,  and 
linearly  polarized  SH  motion  from  the  Love  waves  will  be 
present  in  the  horizontal  plane.  As  the  SH  polarized  Love 
groups  exhibit  higher  rectilinearity  than  the  Rayleigh 
waves,  the  effect  of  a  particle  motion  filter  which  enhances 
motion  displaying  high  rectilinearity  will  be  to  preserve 
the  SH  phases,  but  attenuate  the  vertical-radial  Rayleigh 
wave  motion.  This  effect  is  seen  in  the  filtered  seismo¬ 
grams  of  figure  3.7.  Since  the  surface  wave  trains  are 
dispersive  however,  different  group  maxima  will  arrive  at 
different  times  so  that  the  spectral  content  of  each 
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component  of  ground  motion  over  a  particular  time  interval 
might  be  used  to  separate  Love  and  Rayleigh  wave-trains. 

Simons  (1968) ,  describes  a  process  in  which 
frequency  filtering  using  measurements  on  particle  motion  to 
shape  the  filter  response  is  used  to  enhance  the  signal  to 
noise  ratio  of  surface  wave  groups.  The  discrete  Fourier 
series  of  the  three  components  of  ground  motion  are  computed, 
and  then  each  harmonic  is  treated  independently.  Amplitude 
coefficients  at  each  frequency  are  weighted  according  to  how 
closely  the  three  dimensional  particle  motion  trajectory  at 
that  frequency  corresponds  to  theoretical  Love  or  Rayleigh 
wave  patterns  arriving  from  a  specified  direction.  These 
modified  coefficients  are  then  used  to  reconstruct  the  trace 
in  the  time  domain.  In  the  band  around  the  microseismic 
peak  (5-8  seconds  period) ,  interfering  effects  of  the  back¬ 
ground  will  cause  attenuation  at  those  frequencies,  but  the 
remainder  of  the  spectrum  will  be  enhanced  by  the  filter 
process.  Band-pass  filtering  to  eliminate  the  microseisms 
can  of  course  be  applied  before  the  Fourier  transforms  are 
computed . 

4.2  The  Filter  Function 

If  the  data  traces  are  rotated  as  in  equations  3.1 
and  then  considered  over  simultaneous  time  segments  of 
length  NAt,  where  At  is  the  sampling  interval,  then  for 
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each  component  of  ground  motion,  the  amplitude  and  phase  of 
each  harmonic  is  given  in  terms  of  the  discrete  Fourier 
coefficients  a(nf)  and  b(nf)  by 


Ai(nf)  =  [ai(nf)  +  b^Cnf)] 


b^  (nf ) 

4(nf)  =  arctan  a ; 


(4.1) 


where  i  =  Z,R,T  corresponds  to  the  vertical,  radial  or 
transverse  component  of  motion  and 


n  =  0,1,2, 


N 


Nf 
2 


—  =  Nyquist  frequency  = 


2  •  At 


f  =  4- 


1 

NAt 


(4.2) 


=  fundamental  harmonic  of  the 


Fourier  series 


A  measure  of  the  apparent  horizontal  azimuth  can  be  obtained 
from  the  function 


3  (nf )  =  arctan 


\  (nf ) 
Ar (nf ) 


(4.3) 


and  an  estimate  of  eccentricity  of  the  particle  motion 
ellipse  is  represented  by 


ip  (nf )  =  arctan 


A  (nf ) 
Az  (nf ) 


(4.4) 


=  [AR (nf )  +  At (nf ) ] 


A  (nf ) 
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The  geometrical  relationship  represented  by 
equations  4.3  and  4.4  is  shown  in  figure  4.1.  We  see  from 
the  diagram  that  if  motion  at  a  particular  frequency  is  pure¬ 
ly  radial  then  AT(nf)  =  0  implies  3(nf)  =  0;  while  motion 
which  is  purely  transverse  (A  (nf)  =  0)  gives  3(nf)  = 

K  z 

As  far  as  eccentricity  is  concerned,  if  motion  is  in  the 
horizontal  plane  only,  A  (nf)  =  0  and  ip(nf)  =  while  a 

LS  Z 

predominant  vertical  component  (A(nf)  ^  0)  yields  \p(nf)  =  0. 
If  A(nf)  =  A  (nf) ,  the  eccentricity  exhibits  no  preferred 
direction  and  ip(iif)  =  j.  The  phase  difference  between 
vertical  and  radial  components  is  also  determined  as 

cl  (nf )  =  $R(nf)  -  $z(nf)  (4.5) 

Functions  of  a,  6  and  ip  are  then  used  to  weight  the  ampli¬ 
tude  coefficients  according  to  the  following  equations 

A' (nf)  =  A  (nf) •cosM[3 (nf) ] *cosK[^ (nf)-0] -sinN[a(nf) ] 

u  A 

A'(nf)  =  Ad  (nf )  •  cosM  [  3  (nf )  ]  *  cosK  [i/>  (nf ) -0  ]  •  sinN  [a  (nf )  ]  (4.6) 

R  K 

A^(nf)  =  At  (nf )  *sinM[3  (nf )  ]  *smK[i|j  (nf )  ] 

The  exponents  M,K  and  N  are  determined  empirically.  We  see 
that  the  Z  and  R  components  receive  identical  weights,  and 
that  each  function  of  a,  3  and  ip  lies  in  the  range  (0,1). 


Keeping  in  mind  the  limiting  cases  of  3(nf)  and 
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Figure  4.1  A  representation  of  the  amplitude  coefficients 
calculated  for  the  three  orthogonal  components  of 
ground  motion  illustrating  how  measures  of  apparent 
horizontal  azimuth  and  eccentricity  of  the  particle 
motion  trajectory  can  be  obtained  for  each  frequency 


component . 
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i[i(nf)  discussed  above ,  it  is  evident  that  if  at  a  particular 
frequency  motion  in  the  horizontal  plane  is  purely  radial 
(B(nf)=0),  A  (nf)  and  A  (nf)  are  preserved  and  Am(nf)  is 
attenuated.  Motion  which  is  purely  transverse  in  nature 
(3(nf)=j)  will,  on  the  other  hand,  be  attenuated  in  the 
vertical-radial  plane.  The  parameter  0  in  the  equations 
for  A' (nf )  and  A'(nf)  can  be  chosen  so  that  a  particular 
horizontal/vertical  displacement  ratio  for  the  particle 
motion  trajectory  is  preserved.  If  ijj  ( nf)  =  0  then  the  radial 
and  vertical  amplitudes  receive  unit  weighting  from  the 
function  cos  [ip(nf)-0].  For  motion  which  is  purely  hori- 

n  k 

zontal  ijj  ( nf)^,  and  the  function  sin  [ i/j  ( nf )  ]  applies  unit 
weight  to  the  transverse  amplitude.  The  function 
sinN[a(nf)]  will  attenuate  the  radial  and  vertical  ampli¬ 
tudes  by  an  amount  which  varies  from  1  to  0  according  to 
how  closely  the  phase  difference  between  radial  and  vertical 
components  departs  from  the  theoretical  ^  value  for  funda¬ 
mental  retrograde  Rayleigh  waves. 

The  combined  effect  of  these  weighting  factors  is 
thus  seen  to  enhance  pure  Rayleigh  or  pure  Love  waves  of 
some  particular  period  arriving  at  the  recording  station. 

If  for  example,  motion  in  the  horizontal  plane  is  predom¬ 
inantly  radial,  the  horizontal/vertical  displacement  ratio 
corresponds  to  the  specified  0,  and  the  phase  difference 
between  radial  and  vertical  components  is  j,  then  the  Z 
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and  R  coefficients  will  be  preserved  and  the  transverse 
component  attenuated.  This  corresponds  to  the  case  of  pure 
Rayleigh  motion.  Conversely,  a  dominant  arrival  at  some 
particular  period  on  the  transverse  component  along  with 
little  or  no  amplitude  on  the  vertical  record  corresponds 
to  a  Love  phase  so  that  the  amplitude  coefficient  for  the 
transverse  trace  will  be  preserved  and  the  Z-R  motion  will 
be  attenuated.  Particle  trajectories  which  lie  between 
these  limits  will  be  subject  to  varying  degrees  of 
attenuation.  It  is  evident  that  this  process,  since  it 
considers  the  particle  motion  patterns  over  all  harmonics 
contained  in  the  time  window  chosen  to  compute  the  Fourier 
coefficients,  will  not  be  efficient  if  surface  Rayleigh 
and  Love  waves  of  similar  periods  and  comparable  amplitudes 
arrive  in  the  same  time  segment.  In  this  case,  the  inter¬ 
fering  effects  of  the  simultaneous  arrivals  will  result  in 
a  mutually  confusing  particle  motion  pattern  so  that  little 
or  no  output  from  the  filter  will  be  obtained.  A  further 
condition  on  the  Z-R  amplitudes  can  be  specified  as 

sin  [a  (nf )  ]  =  0,  n<a  (nf )  <211  (4.7) 

so  that  prograde  Rayleigh  type  motion  will  also  be 


attenuated . 
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4.3  Application  of  the  Filter 

A  Fortran  IV  program  has  been  written  which  can  be 
used  to  process  the  three  component  records  from  the  digital 
system  described  in  Chapter  2.  In  practice,  the  Fourier 
coefficients  are  determined  over  the  specified  window  and 
then  modified  using  the  weighting  functions  as  in  equation 
4.6.  As  no  modification  is  applied  to  the  phase  angles  of 
each  Fourier  component,  these  are  used  along  with  the 
modified  amplitude  coefficients  to  reconstruct  the  trace  in 
the  time  domain.  The  window  is  then  moved  down  the  record 
some  fraction  of  the  original  window  length  and  calculations 
are  repeated.  The  final  output  is  the  arithmetic  average  of 
all  overlapping  time  segments  at  any  particular  point  on  the 
seismogram. 

Fourier  transforms  and  their  inverses  were  computed 
using  a  high-speed  algorithm  developed  by  Gentleman  and 
Sande  (1966).  It  was  found  that  a  cosine  taper  applied  to 
about  the  first  and  last  10%  of  the  data  points  in  each 
window  was  required  to  reduce  slight  phasing  effects 
introduced  by  the  discontinuities  in  the  original  time 
series  at  the  beginning  and  end  of  each  window.  With  a 
window  length  of  512  points,  an  increment  of  16  points  (j^) 

was  found  to  give  the  best  results.  As  the  sampling 
interval  of  the  data  was  0.12  seconds,  this  corresponds  to 
a  time  window  of  61.4  seconds  with  a  1.9  second  overlap. 
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The  program  was  arranged  so  that  any  window  from  128  to 
1024  points,  in  multiples  of  2,  as  well  as  the  increment, 
could  be  specified  on  card  input.  The  factors  of  2  were 
required  for  application  of  the  particular  fast  Fourier 
algorithm  that  was  used.  It  is  evident  that  choice  of 
window  length  should  be  consistent  with  the  spectral 
content  of  the  seismogram  over  the  region  of  interest  so 
that  the  fundamental  harmonic  of  the  Fourier  series 
corresponds  to  at  least  a  single  cycle  of  the  longest 
period  present  on  the  record.  Following  the  suggestions  of 
Simons  (1968),  values  of  M,K  and  N  (equations  4.6),  were 
chosen  as  8 ,  8  and  4  respectively,  and  0  corresponded  to  a 
horizontal/vertical  displacement  ratio  of  0.8. 

A  portion  of  a  seismogram  from  a  magnitude  5.4 
earthquake  which  occurred  in  the  Gulf  of  California  on 
March  21,  1969  is  shown  in  figure  4.2  before  and  after  this 
filter  has  been  applied.  The  data  was  digital  band-pass 
filtered  0.007-0.2  Hz  (5-150  seconds)  before  rotation. 
Predicted  arrival  times  of  LQ  and  LR  for  a  surface  focus 
earthquake  are  indicated  on  the  processed  records.  Although 
good  Love-type  arrivals  are  evident  on  the  transverse  com¬ 
ponent,  the  vertical  and  radial  traces  show  little  signal 
except  for  one  well-defined  group.  It  is  thus  apparent  that 
motion  in  the  vertical-radial  plane  did  not  satisfy  the 
theoretical  pattern  chosen  for  surface  Rayleigh  waves,  in 
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Figure  4.2  Portion  of  a  seismogram  processed  with  the 

surface  wave  discrimination  filter.  For  this 
record:  M,K  =  8;  N  =  4;  0  corresponds  to  a 

horizontal/vertical  displacement  ratio  of  0.8. 

A  61.4  second  window  with  a  1.9  second 
increment  was  used.  Digital  band-pass  filtering 
0.007-0.2  Hz  was  applied  to  the  data  before 
rotation.  The  event,  of  magnitude  5.4,  occurred 
in  the  Gulf  of  California  on  March  21,  1969. 

A  =  22° . 


GULF  of  CALIFORNIA  march  2i,i969 


Z  COMPONENT 


l  MIN. 


5hIOm 


-*■' — — ■AAA/W\Ah — 


RADIAL 


RADIAL 


LR 


— ■'W\a/\AAAA/W\Ma' — ./n/Wv^ - 


LQ 


70 


particular,  the  expected  horizontal/vertical  displacement 
ratio  of  0.8.  An  examination  of  the  unfiltered  Z  and  radial 
components  suggests  that  the  apparent  ratio  might  be  >1. 
Interpretation  of  this  observation  has  two  possibilities: 

1.  Few  well-defined  Rayleigh  wave  groups  are 
arriving  from  this  particular  event. 

2.  Conditions  near  the  recording  station,  or 
the  average  effect  over  the  entire,  or  at  least 

a  large  portion  of  the  crustal  waveguide,  are  such 
that  the  major  axis  of  the  particle  motion  ellipse 
for  the  Rayleigh  wave  groups  is  inclined  at  some 
large  angle  to  the  vertical.  This  would  result 
in  an  apparent  horizontal/vertical  displacement 
ratio  quite  different  from  the  specified  value 
of  0.8. 

A  more  definite  interpretation  could  of  course  be  obtained 
by  comparison  of  seismograms  taken  with  an  array  of  record¬ 
ing  instruments. 

The  example  however,  does  illustrate  that  this 
filter  can  be  successfully  applied  to  the  longer  period 
surface  wave  groups.  Data  from  an  array,  once  processed  by 
this  method,  might  then  be  digitally  filtered  with  narrow 
pass-band  zero  phase  shift  filters  so  that  group  and  phase 
velocity  studies  of  surface  Love  and  Rayleigh  waves  could 
be  pursued. 
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Conclusion 

In  the  course  of  this  study,  it  has  been  shown  how 
digital  techniques  can  be  applied  to  three-component  seis¬ 
mograms  for  enhancement  of  the  signal  to  noise  ratio.  If 
the  spectral  characteristics  of  signal  and  noise  are 
different,  zero  phase  shift  band-pass  filters  can  be  effect¬ 
ively  used  to  improve  the  quality  of  the  records.  If 
signal  and  noise  exhibit  similar  spectral  characteristics, 
the  availability  of  three  component  seismograms  permits 
determination  of  filter  functions  which  use  the  polarization 
properties  of  both  body  waves  and  longer  period  surface 
waves  to  enhance  the  arrival  of  various  phases.  Two  such 
filters,  one  particularly  suited  for  the  study  of  body 
waves,  and  the  other  which  can  be  applied  as  a  surface  wave 
discrimination  process  have  been  considered  in  Chapters  3 
and  4.  With  three  component  data  recorded  directly  onto 
digital  tape  using  broad-band  instrumentation  as  described 
in  Chapter  2,  seismic  energy  covering  a  wide  spectrum  can 
be  obtained  for  immediate  application  of  digital  processing 
techniques.  It  is  hoped  that  future  applications  of  the 
polarization  filters  studied  to  data  taken  with  an  array  of 
these  instruments  will  prove  valuable  in  the  interpretation 


of  seismic  events. 
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A.  1 


APPENDIX  1 


A1 . 1  The  Amplifier  Transfer  Function 


A  circuit  diagram  of  the  analog  amplifier  described 


in  Chapter  2  is  shown  in  figure  A1 . 1 .  A  derivation  of  the 
complete  transfer  function  is  presented  by  considering  the 
separate  transfer  characteristics  of  each  operational 
amplifier  stage.  The  input  resistance  of  the  first  amp¬ 
lifier  is  determined  from  technical  specifications  for  the 
Willmore  Mark  II  seismometer  as  given  in  the  Hilger  and 
Watts  operational  manual  (1964).  If  n  is  some  multiple  of 
the  coil  resistance  of  the  seismometer  then  we  have  for  the 
Willmore  Mark  II 


(Al.l) 


where  T  is  the  natural  period  of  the  seismometer  and  D  the 
damping  factor.  With  T  =  1.5  seconds  and  D  =  0.6,  n  =  3.25. 
With  a  coil  resistance  of  3300fi,  the  input  resistance  for 
the  first  stage  is  given  by 


R. 


3300  x  3.25 


m 


(A1.2) 


10.72  Kft 


so  that  R. 


10  is  used. 


in 
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Figure  A1 . 1  Circuit  diagram  of  the  analog  amplifier 

for  the  Tripartite  system. 


' 


10nf 
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For  simplification  of  the  development  that  follows, 
use  of  an  operational  amplifier  with  differential  input  as 
an  inverting  amplifier  will  be  represented  as  follows 
(Burr-Brown,  1963) 


Figure  A1 . 2  General  Form  of  the  Inverting  Amplifier 

with  transfer  function 


(A1.2) 


Zjr  and  Z;  are,  in  general,  complex  impedances.  As  a  matter 
of  notation,  and  when  necessary,  output  from  a  particular 
amplifier  stage  i,  will  be  represented  by  Eq^,  and  input  by 


E-.,  E , . . . .  The  same  convention  will,  hold  for  circuit 
U  2 1 

components  etc. 
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Amplifier  1 

Following  figure  A1 . 2  and  equation  A1 . 2 ,  we 
have  for  the  first  stage 


zf  = 


R  /it oc  , 
eg  J  ' 

R  + 
eq  g  coc , 


so  that 


R 


eq 


1  +  jtOT 


Z.  =  R. 

1  m 


o  i 


i  i 


R 


eq 


R. 


m 


1  +  joox 


T,  =  R  C 
i  eq  i 


(A1.3) 


(A1.4) 


R  is  determined  from  considerations  of  D-C  feedback  in 
eq 

the  amplifier  stage  with  the  following  configuration 
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for  which 


in 


R  2  +  R  3 
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3  J 


E 


o 


eg 

in 


(A1.5) 


Substituting  the  values  of  the  actual  circuit  components  as 
given  in  figure  Al.l  we  have  by  A1 . 4  and  A1 . 5 


R_  _  =  23.1 
eq 

2 

t  =  1.85  x  10”  seconds 


R.  =  10  Kg 
m 

EQi  -  2.31  x  10 3 


Elx  1  +  1.85  x  10  2  jw 


(A1.6) 


The  transfer  function  has  a  high  frequency  cut-off  at  m  = 
8.6  Hz  . 


Amplifier  2 

The  second  stage  is  a  low  pass  active  R-C  network 
filter  with  a  corner  at  0.033  Hz  (30  second  period)  and  an 
attenuation  rate  of  12  db/octave.  Following  the  design 
criterion  described  by  Sallen  and  Key  (1955) ,  the  transfer 
function  of  this  stage  is  of  the  form 


-H|  H 
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E12  ( j <jo )  2  +  d  ( jui)+  1 

K  =  gain  factor 

which  for  this  stage  is  given  by 


l  2  i  +  - 

20.9x10""  20.9x10  2 


(A1.7) 


(A1.8) 


Amplifiers  3  and  4 

Two  cases  A  and  B,  indicated  by  the  circuit  diagram, 
will  be  considered. 


Case  B 

We  have  at  stage  3,  two  inputs,  so  that  the  amplifier 
is  used  as  a  voltage  adder.  Following  equation  A1.2,  the 
transfer  function  for  this  stage  can  be  written 
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03 
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E , 
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1+  joOT 


(A1.9) 


If  T  ( j a) )  represents  the  right-hand  side  of  equation  A1.8 
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then  E  =  T  ( j  go )  -  E  so  that  A1 . 9  becomes 

2  3  2  1  2 


-  R 


f  3 


O  3 


1  + jtOT 
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1  3 
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where  we  have  used  E  =  E  =  E  (output  of  first 

1  3  12  0 1 

At  amplifier  4  we  have 


04 


1  4 


R 


2  4 
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,  T  =  R_  C 

4  f  4  4 


Combining  A1.10  and  Al.ll  with  E  =  E  we  write 
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Substituting  values  of  circuit  components 


R  =  226KQ,  t  =  t4=  10  seconds 


Rf  3  =  lOOKfi  R23  =  R24  =  Rf4  =  lOKfi 


equation  Al.12  yields 


0.443 


10 


+ 


1+10“ 4  (joo)  1+10  (joo) 


T.,  (jm) 


1+10  4  (joo) 


(A1.10) 
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Case  A 


Following  a  similar  development  as  for  case  B,  the 
circuit  configuration  at  amplifier  3  now  gives 
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Amplifier  4  is  now  in  a  summation  mode  so  that 
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Substituting  the  values  of  circuit  components: 


R, ,  =  22.6Kfi  t  =  t  =10  seconds 

14  3  4 


R.  =  lOOKft 

f  3 


R,  =  R  =  R.  =  lOKft 

2  3  2  4  f  4 


equation  A1.15  yields 


-  0.443 


10 


T,  (joo) 


1  +  10  4  (  jco)  1+10  (  j  03  )  1+10  4  (  j  03  ) 


Eol  (A1.16) 


4 

Contributions  from  the  term  10~  ( j oj )  which  corresponds  to  a 
high  frequency  corner  at  about  1500  Hz  are  negligible  so 
that  equations  Al.13  and  A1.16  may  be  written 
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•Ei,  Case  B 
o 


(Al. 17) 


E0<+  -  £-0.443  +  10  T  2  ( j  w)  J  •  E  i  ,  Case  A 


We  see  then  that  case  A  corresponds  to  "out  of  phase" 
summation  and  case  B  to  "in  phase"  summation  which  can  be 
used  for  microseismic  attenuation. 

Amplifiers  5  and  6 

Amplifiers  5  and  6  are  aliasing  filters  each  with  a 
double  corner  at  4.3  Hz.  Their  transfer  functions  are  the 
same  type  as  that  given  in  equation  A1.7  so  that  for  these 
two  stages  we  have 


2 


1.7 


(Al. 18) 


,  ,  1 . 2  ,  .  .  ,  j  oo 

1  +  — 2 ~  (jw)  t  2 ~y 


Combining  equations  A1.6,  A1.13,  A1.18  and  the  expression  for 
T 2  ( j oii )  we  have  the  complete  transfer  function  for  case  B. 
Equations  Al.6,  A1.16  and  A1.18  yield  the  complete  transfer 
function  for  case  A.  These  expressions  are  given  in  figure 
Al.3.  An  APL  program  has  been  used  to  evaluate  the 
amplitude  and  phase  responses  of  the  complete  transfer  functions 
and  these  are  shown  in  figure  A1.4. 
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Figure  A1.3  The  complete  transfer  functions  for  the 

circuit  diagram  of  figure  Al.l.  Case  A 
is  the  case  of  "out  of  phase"  summation, 
Case  B  the  case  of  "in  phase"  summation 
for  microseismic  attenuation. 
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Figure  A1 . 4  Calculated  phase  and  amplitude  responses 

for  the  transfer  functions  shown  in  figure 
Al . 3 .  The  two  cases,  A  and  B,  are  indicated. 
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A1.2  Theoretical  Asymptotes  for  Seismometer  and  Amplifier. 


Consider  the  following  situation: 


Let  x  be  the  position  of  the  center  of  mass 
of  a  seismometer  relative  to  an  inertial 
reference . 

Let  y  be  the  position  of  a  point  on  seis¬ 
mometer  frame  relative  to  the  same  inertial 
reference . 

K  =  spring  constant 
M  =  mass  of  seismometer 

B  =  sum  of  viscous  and  electrical  damping 


factors 


The  second  order  linear  differential  equation  for  the  motion 


of  the  seismometer  mass  is  given  by  (Richter,  Chapter  15) 


1  d2x  2C  dx 


oo  2  dt2  oo  dt 


+  — 


+  x 


(A1.19) 


oo  dt 


K  dt2 


n 


n 


seismometer 


B 


2/MK 


Adopting  Laplace  transform  notation  (Dettman,  Chapter  6) 
so  that  L  [ f  (x)  ]  =  F  ( s )  ,  s  =  joo,  we  have  from  equation 


A1.19 
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s2  2  C  M 

-  X (s )  +  —  sX(s)  +  X (s )  =  -  s  2 Y  (s )  (A1.20) 

to  2  OJ  K 

n  n 


The  output  voltage  generated  by  the  seismometer  coil  is 
given  by 

dx 

e  =  -  G  —  (A1.21) 

dt 


where  G  is  the  transducer  constant.  In  terms  of  the  Laplace 
transform,  this  can  be  written 


E  (s)  =  -  G  ( s )  sX  (s ) 


(Al .22) 


Combining  A1.20  and  A1.22  we  have  as  the  transfer  function 
of  the  seismometer: 
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Considering  velocity  sensitivity,  we  are  interested  in 
Eq(s)/s  Y (s )  so  that  A1.23  yields,  with  s  =  jw. 
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This  expression  may  be  easily  evaluated  for  the  three 
cases  oj < < oo  ,  u  =  a)  ,  a 3 > > oo  to  obtain  the  theoretical  asymp- 
totes  for  the  velocity  sensitivity  of  the  seismometer.  We 
have 


T  ( joo) 
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GM  , 
—  co 
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GM 
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aJ<<0°n  (l^db/octave  slope) 

co  =  oo^  (a  point)  (A1.25) 


T(jco) 


>  G 


oo>>co  (a  constant) 
n 


Considerations  of  the  theoretical  asymptotes  for  the 
amplifier  show  that  the  two  aliasing  filters  give  a  cut¬ 
off  of  80db/decade  from  a  4.3  Hz  corner  and  the  filter  in 
stage  2  has  a  double  corner  at  0.033  Hz  (40db/decade  slope). 
Figure  A1.5  shows  the  asymptotic  response  of  the  seismo¬ 
meter,  amplifier  and  the  combination  of  the  two.  The 
figure  serves  only  to  illustrate  the  bandwidth  available 
from  this  system  so  that  the  amplitude  scale  is  not 


absolute . 
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Figure  A1 . 5  Theoretical  asymptotes  for  the  velocity 

sensitivity  of  the  seismometer  and  the  trans¬ 
fer  function  of  the  amplifier  showing  how 
the  two  curves  are  combined  to  attain  the 
desired  bandwidth.  S  marks  the  curve  for 
the  seismometer,  A  the  amplifier,  and  S-A 
the  seismometer-amplifier  combination. 
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i  1 

2  6 

) 

UR 

I  fL 

(  6 

7 

f  i- 

26 

) 

ST 

UP 

:ND 


SUL-P. 

OUT  I 

NE  L 

r\ 

DEL 

S( 

>11 

N, 

TP  I  ) 

LABEL  P 

LUTS 

COMM 

UN/O 

AT r/  S 

HUM 

TM 

,  0 

AY 

,  YEAR 

,  T 

I  H 

E 

CALL 

S  Y  i : 

C'  Ol  ( 

L 

•  c  * 

3  . 

r; 

■>  f 

^  • 

3  C  ,  '  T 

RA 

c 

!  .  O 

VER 

S  t 

:  •  ,o.o 

,  1  0 ) 

CALL 

SYM 

LOL  ( 

l 

» o  * 

6  . 

v.  T 

n 

L.  9 

3  C  ,  '  T 

A  D 

I  A 

L '  , 

o 

.0,6) 

CALL 

SY  f- 

B  CL  ( 

c 

•  o , 

8 . 

5  } 

w  • 

3C  ,  '  M 

CD 

UE 

ATE 

D 

Z  '  ,  c . 

0,11) 

CALL  SYMBOL  (C.L  ,  11  .C,C.3C  ,  'TRANSVERSE  DIREST  I  (IN' 

0  ,13.9,0.30,'  RAD  I  AL  0  1  RECT  i  DM  '  ,0.0 
0,16.0, 0 . 3 C  t  '  2.  t)  1  RECT  I  ON '  ,0.3,11) 
.C , 18  .  b  ,G .33  ,  'RECT  I  L  I  NEAR  I  TY*  ,C  .0  , 


C  A L  l  SYMBOL  (  C 
CALL  SYMSTUC 


:  v  r‘p 


:l 


C  A  L  L 
C  A  L  L 
CALL 
CALL 
C  A  I.  L 
CALL 
CALL 
C  A  L  L 
C  A  L  L 
CALL 
CALL 


CALL 
CALL 
CALL 
CALL 
CALL  SY 
RETURN 
END 

SUBROUTINE  1  EADSM  LREC  ,  KKK  ,  SC  AL  ) 

LEADS  OAT  A  FPL!)  DISKS  FOR  CO/. PUT  AT  IONS 
TOT  DLLS  AT  ONE  TIME 

COMMON/ ST  OK AG/DAT ( 3,4096 ) 

COMMON /LLK/ I 8LK 

CuLMiCN/!'  I  S  K  S  /  C  h  A  Ml  (  2  3  4  8  )  ,CHAN2(204b)  ,CHAM3(  2048) 
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;b 
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0 

♦ 
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t 

7 
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C 

F 

j 

• 

T 

J 

F 
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A 
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A 

L 

i 

F 

n 

» 

r 

f 

6) 

N 

O 

y; 

*  n 

1  L J 
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( 

• 

c 

J 

2 

6 

• 

u 

F 

ij 

• 

9 

v_ 

F 

•  Z 

c 

r\ 

E‘i 

P 

■m 

H  N 

T 

1 

r  0 

• 

0 

F 

1 1  ) 

5 

V 

'  JO 

1  J 

OL 

{ 

A 
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• 

F 

p 

7 

• 

Q 

f 

■> 

• 

5 

F 
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K 

I 

Q 

i 

A 

R 
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t 

F 

0  . 

j 

* 

1 

0) 

I'i 

1  jl. 

V  ' 

:u 

ER 

( 

• 

F 

2 

7 

• 

r 

f 

pt 

• 

3 

F 

S  M 

ru 

:vJ 

N  ) 

Til 

, 

0  • 

F 

- 

I 

) 

c 

V'  * 

4  ! 

,  8 

OL 

/ 

\ 

3 

• 

b 

’j 

F 
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7 

• 

C 

r 

• 
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F 

97 
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rv 
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-1 ) 
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1  M 

Uf- 
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i  i  3 

ER 

< 

lo 

• 

8 

L. 

F 

2 

7 

• 

I  > 

t 

o 

• 

9 

-> 

u 

F 

DA 

Y 

F 

'J 

• 

0 

F 

1  ) 

c 

YE 

!C 

OL 

( 

o 

• 

'X 
_ * 

V- 

? 

/_ 

7 

» 

0 

0 

• 

4 

'J 

,9 

7 

F 

c. 

• 

c 

t  ~ 

1  ) 

f  « 

Ju 

Uf 

.  B 

LK 

( 

• 

/ 

C 

f 

jU. 

7 

• 

> 

• 

9 

P. 

'•y 

» Y 

r 

AR 

F 

o 

r 
•  Cr 

9  “ 

i 

JL 

) 

•J 

Yf 

■  2 

QL  ( 

7 

« 

1 

0 

F 

2 

~7 

f 

• 

y 

c 

• 

4 

F 

97 

F 

r 

v> 

m 

0 

» 

-1  ) 

u 

Uf 

m 

EP 

( 

7 

( 

• 

z, 

/> 

0 

t 

y 

7 

• 

c 

y 

c 

• 

9 

ij 

fTI 

M 

E 

F 

p; 

r 

•  c 

,  ~ 

1 

) 

SYf 

:B 

GL 

{ 

O 

• 

4 

C 

F 

’  / 

7 

• 

u 

y 

c 

• 

7 

p 

• 

F 

G 

M 

i 

f 

, 

L.  • 

0  , 

s* 

-> 

) 

f  _ 

IJ  f 

.CLP 

( 

9 

• 

4 

c. 

t 

£_ 

7 

• 

9 

0 

• 

r 

F 

TP 

i 

F 

o 

• 

0 

,  2 

) 

\  ! 

10 

OL 

( 

1 

V/ 

• 

Q 

t 

7 

• 

9 

c 

F 

r 

• 

O  ' 
L-  'v> 

F 

f 

c  r 

SE 

c  /  I 

NCH 

• 

rs 

y  U 

• 

0 

F 

8  ) 

N 

Uf 

IB 

E  R 

( 

9 

• 

4 

F 

7 

6 

• 

8 

» 

0 

• 

2  0 

,u 

IN 

F 

3 

• 

■> 

y 

2  ) 

L  (10. 5C  ,26.80,''  .20,  '  SCC  W  INDOW'  ,0.0,  1 


i  v 


KK.FC.OGU  TO  4 


i  I*  l 
L  J  =  K  K  K  -  1 
DU  3  J  J  =  1  ,  L  J 
READ ( 1 ) 

R  L  A  0(2) 

READ ( 3 ) 


,0.0,20) 

,16) 

14) 


2  3 


CUNT  IMUE 
CONTINUE 
1=0 

DO  1  L= 1 , LR  EC 
R l  A ;.)  {  1  )  CHAN1 
RE AO (2)  CHAN2 
R E  A D ( 3 )  CHAN 3 
DU  2  J=1 ,  I  BLK 
1  =  1+1 

DAT (1,1 )  =  C 1 :  A  N 1 ( J )/SCAL 
GAT (2,1) =CH AN2 ( J ) /SCAL 
DAT ( 3, I ) =CN AN3 ( J)/SCAL 


2  CG 

u  r  i 

U  E 

1  cu 

NT  I 

b  » 

UE 

RE 

i  T  M 
» .  1.4 

D 

1 

RE 

W  I  N 

0 

2 

>  r 

WIN 

S J 

3 

."!*  P 

TUR 

N 

.  i_  i 


END 

SUBROUTINE 
C  I  F 


ovrs 


REMOVE ( LREC) 

DATA  NUT  FILTERED 
COMMON/ ST Hk AG/DAT {  3,409c  J /GPSE  T/GC( 3) 
C G 0 i ; U N  /  C  L K / I BLK 
NTCT=LRCC*I BLK 
DU  1  1=1, N TOT 

G  A  T ( 1,1) =DAT (  1  ,  I  ) - DC (  1 ) 

DAT (2,  I ) = ! ) A T ( 2 ,  I  ) - C C ( 2 ) 

DAT ( 3 , I  ) =  DAT ( 3 ,  I  ) - DC ( 3) 

CONTINUE 
RE  TURN 
•END 

SUB  ROUT  I  Nl 


G A I N  ST ( I  ST  ART , LW IN, I  NO LX ,  IB,  I C ) 


I  E  L 


ZH 


n  c 


CuMNOK/uAI NS/RL ( 
CUMDUI./O LK  /  I  BLK 


IF {  I  STAR  T 


AND  SET  INDICES  FOR 
j  5  0 )  , DSV (2050 

TO  51 


FIRST  PASS 
)  ,  DSH  (  2  G  5  0  )  ,DP(  2050 


N  G  «  C  )  GO 


70 


I  A=  (  E  W  I N  -  1  )  / . 
E  U  !  0  I  =  1,1  A 
RL ( I )=C .C 
0 S V (  I ) =0.0 
US  IK  I )=C .0 
DP ( I ) =  t  . C 
CONTINUE 
I  NU  X  —  0 
I  i=  (  L'.j  IN  +  1  )  /, 
GU  TO  52 
CUNT  INUE 
I  \=I START -1 


IPO 


71  1  =  1  ,  I  A 


RL (  1 )=C  .0 


2  4 


C 


C 


c 


3 S V (  I  ) =0.0 
OStil  I  )  =0 . 0 
OP  (  I  )=C-  .0 
71  CONTINUE 


IfiOLX={  1  START-I  LWIN+1  )/2) 
13= I  START 
5?  am  tinue 
I C  =  I  B  L  K 
RETURN 


END 

SUBROUTINE  SET(Ia 


i-  i 


GTS 


VOICES  F  Ci  R  B  F  U I 


IC  , ININ) 

NNING  UE  CLOCK  CN  SUBSEQUENT  PASSES 


10=1 

IC= (LU  IN-  1  ) /2 

RETURN 


•JD 


S  J  0 R U UTINE  k I N 0 0 W (  INDEX, I  STOP) 

FILLS  VECTOR  OVER  WINDOW  FOR  CALCULATIONS 
CON  HCT./STOR  AG/D  ATI  3,4096  ) 

C OR r'.C 0 /LI  SKS/CHAN1  {  2043  )  ,CHAN2(  2043  )  ,CHAN3(  2' 
CON  UJN/ANCLE/TiiETA 

J  J  =  R 


1 48  ) 


OU  1  J=  INDEX,  I  STOP 
J  J  =  J  J  + 1 

CU'V  1  (  J  J  )  =  C  A  T  {  1  ,  J  ) 

Cii  AN  2  (  J  J  )  =  D  A 1  02, J)  *3  I !';  I  THCT  A)  -DAT  I  3  ,  J  )  -COS  I  THETA) 
CHAN  3  (  j  J  )  =1  A  r  I  3  ,  J  )  *5  I  N(  TIILT  A)  +UAT  (2  ,  J  )  *COS(  THETA) 
1  CUNFIKUE 
RETURN 
END 


SUDRQUT I NE 

c 


CCVAR (L WIN, QUAD ) 


C A L C U L A TtS  C  0 V A R I  AN C E  U A T R 1 X  OVER  SPEC  I F I E D  W I MDO W 
CUM  iOf'./D  I  SK  S/CHAN  1  (  2040  )  ,  CMAN2  (  2043  )  ,  CHAN3  (  2048  ) 


PEAL-4  NNX, MNY, MNZ 
DIN ENSIGN  Q UADI  6 ) 
PTS= FLOAT ( LW IN) 

NNX  =  0 . 0 
M!  1Y  =0  .C 


MNZ =0.0 
DU  1  J  = 1 , L  W  I  N 
N  iX=  NNX  +  C !  1AN3  I  J  ) 
NNY=NNY+CHAN2 ( J ) 
iUNZ  =  Mf:Z  +  CHANl  (  J  ) 


I  CONTINUE 
NNX=MNX/PT  S 
N.  j  Y=MNY/  P  T  S 
•1NZ  =  MNZ/PTS 
VAR  1  =  0  .  C 


VAR 

VAR 


C 


2  5 


r 


C  i'  1  \  I  1  ?  —  p.  /*\ 
1/  '  <  v  it-“5  •  _/ 

CUV  1  3  =  C  .  C 
C  JV  23  =  1  .C 


2 


r.j  2  j=i,Lwir: 

V  4  Ul  =  V  A  R  i-MCH  A  N  3  I  J  )  - !  A  N  X  }  *  (  C I  i  A  N  3  {  J  )  -  M  N  X  ) 

V  A  i'.2-  V  A  R  2  +  (  C  t  i  A  M  2  (  J  )  -MNY  )  YICHAN2  (  J  )  -MNY  ) 

V A R  3  =  V  Ak 3  +  {  C  H  YN  1  (  J  )  - MM  L  )  *  (  C  H AN  1  {  J  )  - MN Z  ) 
CUV 12=CGV 12+ ( CHANG ( J ) -MMX ) * ( CHAM 2 ( J ) -MX 
CCJV13  =  CL-V13  +  (  CH  AN 3  (  J  )  -MNX  )  *  (  CHAN  1  (  J  )  -K 
C  0 V  2  3  =  C R  V  2  3  + ( CHAM 2 I J ) -MNY) * ( CHAM  1 ( J ) -  MM  L  ) 
CO NT  I HUE 


*  1  7 
J  t_ 


ART  j=FI.  LA  T  (  LWIN  ) 

VAR  I  =  VAK1/PPTS 

V  A  'A  2  =  V  A  K  2  /  P  P  T  S 

V  AM  3= VAR  3/PPTS 
CUV  12  =  CUV  12/PPT  S 
C  DV  i  3  =  CCV1  3/PPTS 
C  U  V  2  3  =  C  C  V  2  3  /  P  P  T  S 
uUADI I ) =V Ak 1 
QUAD  (2  J  =  C0V1  2 

QU AU ( 3 ) =VAR2 
QUAD (A ) =  C  0  V 1 3 
QUA!;  (  5  )  =CC)V2  3 
.UADI o) -VAR 3 
RETURN 


S  J 3P  GUT 1 N  2  G A  I N A  V  <  NW I N , X X ) 

AY Ck AGES  THE  GAIN  FUNCTIONS  OVER  A  SPECIFIED 
COM;  UM/3LK/TRLK 

CG'-MCM/CAIf’S/RL  (  2050  )  ,DSV(  2050  )  ,DSH{  2050 
01  MENS  I  ON  XX ( 1C ) 

K=2*IELK 

KK  =  3*  I BLK 

P T$  =  FLOAT ( MW  IN ) 


WINDOW  LENGTH 
DPI  20  50 ) 


I  NO  E  X  =  C 

i\ST=(MWIN  +  l  )/2 
RENU  = I  CL  K-( NWIN-1 ) /2 
MST't  1  =  NST  -  1 
NPLUS  =  NENU  +  I 
DU  1  i  -  N S  T  »  N  E N D 
I^INCEX  +  NWIN 
INDEX- I r  Utx  +  l 
AV1-0.C 


A 

M 


V  2 


C 


A  \r~> 

/  \  i  c 


,nV  .  U 

uN  2  J=IMDtX,IC 
A  V  l  =  A  V 1  +  R  L  I  J  ) 

AY 2  =  AV2  +0  S V ( J  ) 
AV3-AV3  +  D5HI J  ) 


<1 


c 

f 


c 


c 


c 

c 


A  V  4  =  A  V  4  +  D  P  (  J  } 

COMT  I r. U E: 

A  V 1  =  A  V 1  /  P  T  S 
A  V  2  •-  A  V  2  /  P  T  S 
AV 3= AV3/PT  S 
AV4-AVA/PTS 
XX  (  ! ) = A V 1 
XX  (  I  +  i  L  L  K  )  =  A  V  2 
X  X  (  I  +  R )  =  A  V  3 
XX (  I +  K  K )  =  A V 4 
1  CONTINUE 

ill.  GAINS  WITH  AVERAGES 
EG  INN  IMG 


Du 


1=1 i NSTMl 


RL (  I  ) =XX( MST ) 

03V  (  I  )=XX(  NST-HBLK  ) 
OSH{  I )  -  XX  (  MST  H<  J 
OP (  I ) -  XX ( MST+KK ) 

3  CONTINUE 
0 

03  4  I  =  NPLUS , I BLK 
■'.L  (  I  )=XX(  MEND) 

0  3 V  {  I  )=XX  (  NEND+- 1  BLK  ) 
D S i  I  {  I  )  =  XX ( N  END  HO 


J  N!  ' 


rw 


(  I  )  =  X X  (  NEMD  t-KK  ) 


4  CONTI  r :  U  E 
IDOL  E 


DO  5  1  = 

r 

c 

vJ 

T , N  END 

FL( I ) =X 

V 

t  \ 

( 

I  ) 

o  s  v  £  n  = 

X 

V 
/  . 

( I f I BLK) 

DS  II  I  )  = 

X 

y 

(  I+K) 

DP ( I ) =X 

V 

c 

I  +KK ) 

co  pit  i  r ;  u 

c 

RETURN 

E  JD 

SUE  ROUT 

I 

• 

* 

E  ROTATE  ( 3 

.  T,Z  ) 

ATES  GR 

i 

r 
' } 

INAL  DATA  T 

RACES 

,Z  CG-U 

i  V 

D 

IRATE  SYSTE 

M 

INTO  AN 


C  JNULR/3T0R AG/DAT  {  3,4006) 

COOPQN /ILK/ I  BLK  /  AM  jL  6/ THETA 
C  I  'IE  NS  I  L  N  R  (  10)  ,T(  10)  ,Z(  1C) 

33  1  L= 1,1 BLK 

R ( L ) = C A T ( 3 , L ) *S INI  THETA) +  DAT ( 3 , L) *COS 


AT  I 


PPL 


T  (  L  )  =  L  A  7  (  2,  L  )  ;':  S  I  M  (  THETA)  - 
Z { L ) -P  AT (  1 , L  ) 

C  O'  N  T  I  f .  U  E 

RETURN 

Efil) 

SJLROUTINE  GAINtR, T , Z) 

FUMCT I  CMS  TG  THE  DATA 


, L ) NCOS 


-  .3 


O'  /  '  1  * 


( 1 HETA ) 
( THETA  ) 


2  7 


C 

C 

c 

L 


c 

UN 

:n 

./CA 

INS/RL 

(  20 

AN, 

D S  V  (  2 C 5 0  )  ,  D  S  H  (  2  C  5  O  ,DP(  20  50 

c 

')  M 

M  C 

U/DL 

K/  I  OLK 

0 

r  , 

i  i  < 

E  N 

SION 

R  (  1 0 ) 

,71 

10)  , 

Z  (  1  (  ) 

n 

U 

1  I 

L  =  1  , 

I  B  L  K 

Z. 

(  L 

)  = 

Z  {  L ) 

*RL (  L) 

4  D  P 

(  L  ) 

f\ 

(  L 

)=! 

!  {  L) 

*RL  (E) 

*DS 

V (  L) 

y 

l  L 

)  = 

ML) 

*RL( L ) 

*DS 

H(  L  ) 

1  CONTINUE 
RETURN 
END 


r 

J 

UB 

ROUT  IN 

c 

l_ 

READ 

TP  (  N 

3LK  ) 

NEC! 

> 

T 

DATA  T 

PE  UN 

FT  0 

4 

D  I  S 

f\ 

C 

O 

1-3  1-0 

R 

CHAN 

NEWS 

1 

'■> 

-3 

r  LA 

1/ 

c 

■i  CLKS 

D  E  S  I P. 

ED  A 

NO 

PUT 

S  ONTO  DISK 

CAL  CU 

LA 

IE  DC 

r 

01 1 

ION/  T  E 

r/ 

P/  I  J  A 

T  (61 

4  4 

)  /  0  F 

SET/ DC ( 3  ) 

c 

L 

0  i  ’» 

DON /EL  K 

/  I  BEK 

r 

■jr. 

NU1 / 0 I 

KS/Ch 

AMI  ( 

*  • 

C  _/ 

4  8  )  , 

CHAN2( 2048) , C 

DC ( 2)=G.C 
DC { 3 ) =C  .C 

FT  5  =  EL  CAT  {  I  R  L  K  *  ?!  BLN) 
DU  1  L  =  1  ,  NC  LK 
1=0 

R E  A '.) (  A  )  I  DAT 
DU  I  0=1,6144,3 
I  =  H-1 


CriAM  {  I  )  =  FI.  CAT  (  I  DAT  (  J  )  ) 

C  t  i  A  N  2  (  I  )  =  F  L  0  A  T  (  I  OAT  (J  4-1 J  ) 

C H A ‘ : 3  (  I)=FLC  AT  (  IDATt J+2)  J 
DC(  1  )  =DC (  1  ) 4  CHAN  1  {  I ) 

DC  {  2)=LC  (  2)  4-CHAN 2 (  I  ) 

D C  (  3  )  =  L C  {  3)  4-CHAN3  (  I) 

2  CONTINUE 

WRITE (1)  CHAN  1 
WRITE (2)  CHAN2 
WRITE (3)  CHAN 3 

i  continue 

DC  {  i )  =1  C  (  1  )  /  P  T  S 
DC  (  2  )  =FC  (  2  )  / P T S 
DC ( 3  )  =  D C ( 3) /PTS 
WRITE (4, ICC)  (DC  INK) ,NN=1 ,3) 
ICC  FORHAT ( • COf FSETS’ 3E2C . 7) 

FEW  INC  1 
REWIND  2 
REWIND  3 
RETURN 
E  •  J  0 

SUDKOLiTINE  SKIP(N) 

C K I P 6  DOWN  THE  TAPE 


1 


c 

c 


c 


L 


c 


p  r 


DU  1  J  =  1,N 
REA0(4) 

GOUT  I  NUE 
PC  TURN 
END 

SUdROUT INE  CUT ( NDLK , I  Z  $,  XX  ) 

.AOS  DATA  T  RUM  DISK,  3UTTERW0RTH  FILTERS 
RUTS  SACK  UNTO  DISK 
D I  1  ENSIGN  X  X ( 1C) 

C  0  ;  NCN/D  I  S  K  S  /  C  H  AN  1  (  2  0  A  3  )  ,  C  HAN  2!  207  8  )  » Cl  I A  N  3 

CUNNGN/ULK/  I  DLK 

CONf'CT  /OF  SET  /  DC  (3)  /F  I  L/T  1(8) 

I  <$=NCLK#  I  BLK 
IXZ  $  =  1  X  0  +  IZS 
I  I  =  I  X  $  +  1 

r  fi  AN  (  E  ,  80  )  FREQLC  ,  F  l  EQHC  ,  C 
PE  ATM  3,81)  FI 
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x ,  x  r 

,  KKK 

I 

F 

I  UR 

I  T 

* 

o 

-J 

) 

j 

»  .  !  \ 

T 

T 

r- 

( 

t 

»  1 

1 

4) 

WR 

1 

r 

r 

{ 

6 

?  i 

... 

b  )  (XT  (  I  ) 

f 

1 

= 

i  » K 

KK 

) 

w 

¥ 

A 
r \ 

L 

u 

L 

or 

LOG  > 

N  FO 

Q 

\ 

F  F  T ' 

r* 

A 

i  r 

(  L 

?  j 

v . 

\ 

i 

f  ’ 

•  L. 

L 

.  1  2  3  )  L  OG 

2 

Ni 

- 

7 

I  F 

( 

|o 

I 

R 

L". 

•  t. 

,  V 

.256 

)  LOG 

/ 

N 

= 

Q 

IF 

( 

L 

T 

i 

f ; 

.  L 

ft 

• 

Ul 

)  LOG 

K 

i  ^ 

- 

o 

J 

i  F 

( 

L 

, 

I 

f 

i ; 

.h 

? 

}  7 

•  x.  _ 

4  )  L  0 

G 

■p 

N 

—  1  L 

i  XI 

1 1 

= 

1/ 

K 

-L 

U 

IN  +  1 

Cl 

r  • 

L 

■j 

T 

T  H 

ARAL 

YSI5 

II 

r\ 

7 

OR 

N 

I 

NC 

FILL  V 

ZCTCK 

C 

F 

0 

R 

TWO 

KFT 

C  A  L 

L  KIR 

Q 

OK 

( 

T 

F 

lv;  I 

R  ,  XZ 

3"  l' 

A 

fv  * 

c  c 

K  ml1  i  at 

»  ♦  •  f  i 

1 

I  r 

/-i 

1  , 

j 

v  r_ »  * 

ij  (Z  0  • 

FILL  1 

i  < 1A ?  '  V 
j  M  .  1 

r; 

i 

OR 

T 

0  R 

FT 

DG 

1  JJ  = 

1 

>  L 

V 

IN 

XTfXR»T7»FT» FR) 

!  i  L  m  •  »  1  Z 


c 

c 


c 

c 

c 

c 

( 

c. 

c: 


/>.  r> 

t  -\  i 


.CM. 


ccv 
I'D. 
OUT 
TO  A 
I 
i. 


N 

I  I 


DU  i!iY  (  J  J  }  =1 ,0 
CONTINUE 
EK 

T A PEP. ( LW  I N i NO PIS ,00,01,02) 

T  !•  /\ N S 0 OiCI S 
T l.UkFT  {  LuG2N,  0  2  ,  FR) 

01  CO: FT  (  L CG 2 U ,  FT  ,  DU !'!•', V  ) 

PGUP  ILF.  TRANSFORMS  OF 
SEPARATION 

P  U  T  Ff  CM  S  U  00  GUT  I  N  0  S  W  A  YE  IS  M  0  DIP!  £  D  R  £  C  0  \!  STRUCT60 
C  t  !.  •  V  t.  P  W  I  N  D  C  -J  IN  I  I  M  C  U  0  U  AIN  STO  R  6  D 
FK  .  .  .  \ .  A  L)  I  AC  ,  PT  .  .  .  T  KAN  S  VERSE  .02  ...  2 
TE  THESE  SEGN 


CALL 
PU  TE 
C  A  L  L 
C  A  L  L 


!  i  A  V  E 
FORM 


R,T 


!T S  UN  DISK 


C  A  L  L  S  U  A  V  C  (  L  W  I  N  ,  L  U  G  2  N  ,  F  2  ,  F  P. ,  F  T  ,  A  2  ,  A  R  ,  A  T  ,  P  H I  Z  ,  P 1 1 1  P  ,  P  ■ 
a  R I  T  E  (1  )  F  2 
WRITE (2)  FT 
K i .  I T  e ( 3 )  OR 
2  CONTINUE 
I- ii  Vi  INC  I 
REWIND  2 
REWIND  3 

RECONSTRUCT  ENT  IRE  TRACE 

Call  STACK!  NBLKS,L  .IN,  I  IlC  ,  X  Z  ,  XT  ,  XR  ,  F  2  ,  F  T  ,  F  R  ) 
PLOTTING  OUTPUT 
XZ ( KKK  +  2 ) =  AMP  2 
XR ( KKK+2 ) -ANP2 
XT ( KKK+2 ) -AMP 2 
CALL  PLOT ( C .0,-2. 5,-3) 

OKI  T  E ( 6 , 102 ) 

W  P I T C I  6 , 1 C  8 )  ( X Z (  I )  ,1  =  1,  K KK ) 

C  At. L  L  I N  E  I  X  ,  X  2  ,  KKX  ,  1  ,  I  UiN  I  T  ,  3  ) 

CALL  PLOT  I C . C , -2 . 5 , -3 ) 

WhT TE ( 6, 103) 

a;1  I  TE  (  6 , 1  08  )  (  X R  I  I  )  ,1  =  1,  KKK  I 
C  41 L  LINE! X , X  R , KKK , 1 , I  UN  I T , 3 ) 


C  K I  T  0 

(  t 

,  i  c 

4 

) 

U  I 

1 — 

l - 4 

* 

(  6 

1  ~ 
f  I  - 

8 

) 

(XT  (  I 

) , 

r  _ 

i  “ 

1 ,K.KK ) 

CALL 

P  LOT  ( 

• 

»  “  2  • 

5  , 

-3 

) 

CALL 

L  I 

r  E  ( 

X 

> 

XT  ,  KK 

K  J 

1  , 

I  UN  I T  , 

CALL 

PL 

GT  ( 

c 

• 

b  ,  b  •  ~ 

,  9 

99 

) 

WRITE 

(  6 

,  11 

L- 

) 

C  R  I  T  r 

(6 

,11 

2 

) 

STOP 

END 

SUB POL  T  T  N E  STACK ( NELKS ,  L  N  IN  ,  I  NC  , XZ , XT , XR , F2 , FT , FR ) 
STACKS  DAT/.  OCR  PUP. PjsES  OP  PLUTTlNb 
L  1  Mi!  NS  I r  h  X  2  (  1  .  )  ,  X  T  (  10  )  ,  X R  (  1C  ) 
u 1  •  1 E MS  ION  0  2(1024)  , FT (  1  024  ) ,TR(  1024) 

C  0  i  i  H  C  K  /  C  L  K  /  I  BLK 


II  T) 


UP T  o  =  uuLKS*I GLK 
v:  I  r :  C  =  F L  G A  T  (  L  W  I  N  )  / F  L  0  AT  {  I  f'JC ) 
OG  1  1  =  1  t  N  P  T  S 
XZ  l  I  ) =0.0 


X  T  (  I  )  =  C  .  0 
XR (  I) =0.0 
COM  IMUE 

UTJ PTS^ril  TS-LVUN+  l 
1=1 , NNPTS , I UC 


Uu 


E AO ( 1 )  FZ 


R  E  A 

0(2 

)  FT 

F  l  A 

D  (  3 

)  F 

R 

IMD 

ZX  = 

1-1 

• 

DU 

0  J 

=  1  i 

LWI  N 

I  HO 

E  X  = 

IMD 

EX+  1 

x  Z  ( 

i  J  s’  L: 

EX  ) 

=  FZ.  (  J  )  t-XZ  ( 

T 

i 

NDE 

X) 

XT  { 

1 1 .  D 

LX) 

=  F  T  ( 

j )  +  x  r  ( 

I 

NOE 

X) 

xr.  ( 

t 

4.  1  - 

c  X ) 

=  FR ( J ) +  X  R ( 

I 

NDE 

X) 

CUN 

r  n 

UE 

COM 

TIN 

UE 

!_  A 

VET 

AGE 

S 

Du 

4  J 

-1  , 

U  HP  T 

c 

s J 

XZ( 

J)  = 

X  Z  ( 

J  )  /  W 

INC 

XT  ( 

J  )  = 

XT  ( 

J  )  /  W 

INC 

XR  ( 

J)  = 

XI;  ( 

J  )  Z  G 

I M  C 

CUN 

r  i  r- 

UE 

R  E  G 

I  N  D 

1 

REW 

INC 

C 

REW 

INC 

n 

it 

RETUP.f. 

END 

SUP 

RUU 

TIN 

*  1  i 

E  Si. 

AVt l LG 

I 

n,l 

U  0 

P  EPFCRM-S  0  UK  FACE  WAVE  D I  SCR  I  M  I  NAT  I  UK  PRCCCS 
A Fill  1 9c  8  EULL.SE  IS.SOC. AH . 

F  R  ,  r  T  ,  ■  ’ 


...VECTORS  CUNT  A  I N I  NO  FOURIER  TRANS 


A Z  A;  M  .  .  .  CU H P U T t U  A KPLITUDES 

phiz,  run  ,  pH  n...  phase  angles 

CUKHL-f  /EX  PM /IP!,  IK,  IN 
CGKitON/liV/OR 

1  N 


G  F  FT  •  S 


z  ( l  G )  ,  r 


{ ic) *  f  r ( ic) 


( 1C ) , A  (1C), AT { 1C ) 


COM  MOT /PIE/PI , TGOPI 
U  I  HE  NS  I  IN. 

01  KENS  I  Eli 
L)I  Ml  NS  1  L  N  PHI  Z  (  10)  ,  PMIR(  1C 
C  If.  =  FLOAT  (  LOIN) 

NL  P  Jl.c  =  LW  I  r ■  /  2 
r:pLUS  =  f;uPGN2+ 1 

•  L  CUL  AT  L  AMPLITUDE  AND  PHASE 

•  R » "  0  A  IF  FAC  T  C  k  S  AM!)  i  '  f  1 ..)  I  F  Y 


, PHI T (  10  ) 


TERMS 

UjF FT IC  I  ENTS 


r 


re- 


J  —  -  i 
}L  1 


UP1 


? 

L. 


J=2 , NUPON2 


AS,  AT, PH  I Z, PHI K 
S  AS  PER  SIMONS 

FORMS 


PHI  T) 


36 


v.= 


L  ‘  1 i\+2  -  J 


C 

C 


c 


c 

c 

r 

c 

c 


A  L  I  uR  A R  V  S  U  R  F  U N C  T  1 0 N 


O  )  =  1 . 0 


JM 1  =  J  -  1 

AZ  (  J  )  =SQR  T  (  F  Z  {  J  )*  *2+FZ(K)  **2  ) 

AH ( J ) =SQKT { FK ( J ) **2+FR ( K ) **2 ) 

AT  (  J  )  =  £  j  R  T  (  F  T  (  J  )  **  2  f  F  T  {  K  )  **  2  ) 

CHECK  F G  i .  A T  A N 2(0.0, 0 . 0 ) 

Af'.D  FORM  PHASE  AH  CL  L  S  .  .  .  AT  AN?  IS 

I  F  (  {  F  L  (  F. }  .  E 0  .  C  .  C  J  .  AND  .  {  F  Z  l  J  )  .  EG  . 0 . 0 )  ) FZ (  f ’  = 

If(  ( Fi  (K)  .PC  .Z  .  C  )  .  AM D  .  ( F  R { J )  .  L  C  .  0 . 0  )  }  F  R  {  J  )  =  1 . 0 
I F ( <  FT ( K )  . EG.G.C  )  . AMD. ( FT ( J ) .EQ.O.C) )FT( J)  =  I.O 
PHI  l  (  Jf'l  )  =  A  T  A  N  2  (  FZ  (  K  )  ,FZ(  J)  ) 

I  F  (  P  H  I Z  (  J  01  )  .  L  T  .  C  .  j  )  P  H  I  Z  (  J  K  1 )  =  P  H  I  Z  (  J  N  I  )  +  T  WOP  I 
PHIL  ( vif‘1 )  =/,  T  AN2  {  F  R  ( i<  J  ,FR(  J)  ) 

IF [ PH  IF  ( J  T  1  )  ,L f  .C . C ) PHIL ( JM  1) =  PHIR ( JM1 ) +  TWOP I 
PMI T ( J Ml ) = AT A M2 (  FT ( K  ),  FT ( J ) ) 

1  F  (  P  Hi  T  (  J  Ml  )  .  L T  .  0 . C ) P H I T ( J M  1 ) =  P M I T { J M 1 ) +T WOP  I 
:  UL- 1  F  Y  A  OF  L I TUD  E  Cc  EFF  IC  I  ENTS 

2  E T  A- ;  T  AH 2 ( AT ( J  )  , A R ( J  )  ) 

A  =  S  ;  K  T  {  AR  (  J  }  *  *  2  +  A  T  {  J  )  *  *  2  ) 

PS  I  -  AT  AM 2 ( A , AZ ( J  )  ) 

A L F A =  F H I R ( JHI ) -PHI Z  (  J  HI  ) 

IF (  (ALFA. OF. PI)  .AMO. ( ALFA. LE.TWOPI  )  )GG  TO  10 
A 2 ( J ) = AZ ( J  J - (  { CL  5 l MET  A )  ) * *  I M } - (  ( COS ( P  S  I  -OR  )  ) ** IK)* 
.  (  (  S  i  M  A  l  F  A  )  )  **  I  N  ) 

Ail  (  J  )  =  AP  (  J  )  *  (  (  C  0  S  {  B  E  T  A  )  )  *  *  I  M  ) 

.  (  (  Sir.  (ALFA)  )  **1N) 

AT ( J  )  = AT ( J  )  * (  ( S I B { BETA )  ) ** I M ) 

GO  TO  15 
1C  COP.  TIM  UE 
AZ ( J ) - C  .  0 


(  (CCS (PS  I -DR)  )**IK) 
(  (  >  I N  (  P  5  I  )  )  *  *  I  K  ) 


B ( J )=L  .0 


TA)  (SIN  (PS  I)  )  **  IK  ) 


HERMIT  I  AM  FORM 


AT  (  J  )  =  A T  (  J  )  *  (  (SIMS 
15  CONTINUE 
1  COLT  I  HUE 

P LL  0 A D  f !  C LIT  IE 0  AM PLITUDE  COEFFICIENTS  A  NO 
PHASE  ANCLES  INTO  rZ,KP. ,  FT ,  TO  TAKE  INVERSE  FT 
TO  DO  THIS  FORM  REAL  PART  A* COS (PHI )  AND  IMAGINARY 
PART  A- SIT  (PHI)  AMD  STORE  IN  COMP AC 
\S  DESCRIBED  it,  SUBROUTINE  TWO  TFT. 

DU  4  J  =  2  »  NUPQN2 
K= LN  I !  +2- J 
J  M  1  -  J  -  1 

FZ ( J ) =  AZ ( J ) *CuS ( PH  I Z ( JM1 )  ) 
r  Z ( K )  =  - AZ ( J ) * S I R { P  H I Z ( JM 1 ) ) 
n-  (  J  )  -  A F  (  J  )  *CGS  (  PH  1R  (  JM1  )  ) 

FR  {  K )  -  -AF.  (  J  )  *S  I  N  (PHI  R  (  JM1  )  ) 

FT ( J ) = AT ( J ) *  CCS ( PH  IT ( JM1  )  ) 


ASIN(PHI)  IS  COMPLEX  CONJUGATED. 


F  T  (  !'.)=- AT  {  J)*. 


,  (PH  IT  {  JIU  )  ) 


C 


4 

TAKE 


CUNT I  RUE 


INVERSE 


t  .  )>mc  r  i’.D 
i  r>  i  ,i  i  j  Ui\  i  i 


.o 


3  7 


•  CALL  ONEC FT { L0G2N,  FZ) 
CALL  L  iiLCF  T  (  L OG 2 K  »  L R  ) 
CALL  LNECFT ( LGG2N, FT ) 
NOR  '  Ml  If. 

DG  3  J=1 ,  L IN 


T  I  {  j  )  =  f  2  (  J  )  /  W  I N 
LG  U)  =  FR  (  J  )  /  WIN 
r  ri  J )=f  T( J) /WIN 
CON  I INUE 


RETURN 

END 

SUB ROUT  IN 

E  r I LL ( NB LKS ?  XZ  » 

LLS  VECTOR 

S  X Z  t X T j X N  WITH 

TA  TRACES. 

Cl  IAN  1 . . . Z , CHAN 

D  I  IE  NS  I  ON 

X  Z  (  1  j  )  ,  X  T  (  1 0  )  ,  X 

COOOQN/l  I 

SKS/CtiAM  {  2D  48  )  , 

CuN  OON/GLK/  I  BLK 
C O’ILKJf,/ ANGL  E/  THE  T A 
S  =  S  I  r;  (  THETA) 

C  =  CNS (THETA) 

DO  1  J=1 , N3LKS 
RUAO(l)  CHAN  1 
L.  h  .  ( «._  )  C  HA N  c 
KL  \D  (  3  )  CHAN 3 
I  i if.' *.  X-  ICLK*(  J-l  ) 

DO  2  I K= 1,1  ELK 

I  N  D  E  X  -  1 1  i  L)  E  X  + 1 

XZ  (  I ! ;  f.  ■  >. )  =CHAN1  (IK) 

XT  (  If. LEX)  =CHAN2(  IX)  '*$  -C  HAM  3  {  IK)  *C 
XL.  {  INC  F>  )  =  C  H  AN3  (  IK)  *5+CHAN2  (  I  K  )  *C 
2  CONTINUE 
I  CONTINUE 
RE W  INC  1 
RCW  I  NO  2 


REWIND  3 
kE TURN 
END 


c  i  t 
-> 

>) 

i  UUT 

T 

I 

» • 
i  % 

E  N  I 

NO  01  ( 

IS 

TA 

p 

T  ,  LW  I 

NtXZ, 

XT  ,  XF:  ,  FZ  ,  FT,  FR  ) 

F  ILLS 

VECT 

n 

i< 

S  FZ 

,  FT  , 

FR 

i  i 

i 

TO  PC 

I  NTS 

OVER 

W  I N  0  C 

i . » 

PRO 

j 

XZ  r  X 

T ,  X  R 

p  j- 

SR 

CT  I  VE 

LY 

r>  r 

j 

c.  ’.j  r  T 

«* , 

f.  • 

.  i 

X  Z  ( 

Lv  )  ,  < 

T  ( 

10 

) 

,  XL.  (  I 

C  ) 

0  I 

1 

,  ) 

L  N  S  I 

v_ 

4 

FZ  ( 

id)  ,  r- 

T  ( 

l  : 

) 

,  FR(  1 

C  ) 

IN 

{’) 

c  X  =  I 

c 

T 

i 

AP  T  - 

i 

DO 

1  J  = 

L 

7 

LW  IN 

T  *  : 

\j 

EX=  I 

N « 

n 

-  X  +  i 

l_  T 

1  i_ 

( 

J  )  -  X 

2. 

( 

I  Dc 

X  ) 

FT 

( 

J  )  =XT 

{ 

I  NOE 

X  ) 

F  R 

( 

J  )  =  X 

r 

( 

INDL 

X) 

1  CONTINUE 
P.L  TURN 
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C 

r 


C 

C 

c 

c 

c 

c 

c 

c 

r 


C 


c 


EM 

D 

SO 

p 

KCJLT  I 

ME  T 

A 

p  E  K 

ER 

a 

F  IRS 

T  (  I  /N 

FPT 

X  l— 

a 

OF  H 

CO 

I — 

V; 

ITH 

■K 

• 

5  C  0  S  ( 

NFPT 

6 

*P  I 

01 

•  « 
i  i 

ENSIL 

N  X( 

]_ 

0  )  , 

NT 

r  a 

P=RPTS/NF 

P 

TS 

PI 

3  .  1 4 1 

8  92  6 

536 

a  r 

r  \  K 

•  J 

=  FLO  A 

T  (  NF 

P  T  S  ) 

no 

10  JA 

=  1 ,  N 

T 

AP 

JA 

= f :  p  t  s 

—  J  A  + 

l 

FJ 

/  \ 

=  FLOA 

T  (  J  A 

- 

1 ) 

W  = 

0 

•  5  +  C  . 

5  *CQ 

C 

(  A  R 

X( 

J 

A)  =  W* 

X  (  J  A 

) 

Y  (  J 

A  )  =  \-l* 

Y  {  J  A 

) 

Z  ( 

1 

J 

A)=W  * 

Z  (  JA 

) 

X(  JAA)  =  i.*X(  JAA ) 

Y(  JAA)  =  L'*Y{  JAA) 

Z  {  JAA)  =  M*Z  (  JAA) 

LC  CONTINUE 
RETURN 
tMO 

S u G i'  (JUT  IN E  T  V : G  k  F  T  (  LOG 2  N ,  X  ,  Y  ) 


T VIC-  Ri  At.  FOiJl  ICR  TRANSFORMS 

TEL'  EFT  COMPUTES  THE  FFT  OF  TWO  REAL  TIME  SERIES 

X  AN  J  v.  Tilt  our  PUT  IS  STORE!)  IK  COMPACT  HERMIT  IAN  FORM 

THAT  IS,  FOR  N  POINTS  IN  EACH  REAL  TIME  SERIES 

CO  HAVE  R/2  +  1  F\  E  A  L  CO  EFF  IC  I  ENTS  AND  M/2-1 

IMAGINARY  COEFFICIENTS.  FOR  EXAMPLE,  THE  FT  OF  THE  TIME 

X  AS  OUTPUT  FROM  Tub  RFT  KILL  DE  STORED  AS  FOLLOWS: 

X  {  1 )  =  i)C  COEFFICIENT  (  PURELY  REAL) 

X  (M/2  +  1  )=NYLiUI  ST  COEFFICIENT  (PURELY  REAL) 

AND  FOf  J=2,N/2J  K=N+2-J  WE  HAVE 
X  (  J  )  =  k E A L  COEFFICIENT  OF  THE  FT 

X ( K)=CHRr .ESPONOING  IMAGINARY  COEFFICIENT  OF  THE  FT. 
INTEGER  LOG2N 
REAL  X  (1C),  Y  (IN 


INTEGER  J , K , N , N  OVER  2 
REAL  X I , X  R , Y I , Y R 


CALL  HR  10  CT  (  L  CG  2  N  ,  X  ,  Y  ) 

0=2-* LOG 2  N 

N  OVEI  2=0/2+ 1 

DJ  ICC  J  =  2,N  OVER  .2 

K=N+2-J 

XR= Y ( J  )  +X ( K ) 

XI=Y( J)-Y(K) 


SERIES 


C 


39 


YP.=  Y{  J  )  4-Y  (  K  ) 

Y  1  =  X  (  K  }  -  X  (  J  ) 
x(i;j  =  xi#.5 

X ( J  > =  XF  *. 5 

Y  (  K  )  =  Y  I  *  .  t> 

Y  {  J  )  =  Y  K  - .  5 

ilo  cor.  riNuc 
RETURN 
E  f  J  D 

SUBROUTINE  MR.  1  D  r-  J  (  LUG 2 N ,  X  ,  Y  ) 


FIXED  I 


,D 


IX  ONE  LI 


M  C  T  f  >  7  ■ 
V  O  i  '  J  i  M 


L  FOURIER 


t  r  A  FS  P  QRM 


INI EGrI  LUG2K 
R L  A  L  X  (1C),  Y  (10) 
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h=2**LLG2N 

II  (LLG2N.LE.1)  GO  F0  500 
U(0  4CC  K=  2  ,  L  CG  2  N , 2 
0=2**  { LOG  2N-K ) 
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00  30 C  0=1,  M 

AO  G  =  6 . 2  c  3  1  C  5  3  0  7  *  FL  J  A  T  (  J  -  1 )  /  f  L  OAK  H4  ) 
C 1 =C  OS (  A  FIG  ) 

s  i  =  s  i  r  (  a  r<  g  ) 

C  2  =  C 1 *  C  1  -  S  1  *  S  1 
S^  =  C  14SH-C1-S1 
C  3  =  C  2 *  C  1-S2*S1 
S3=C2*S1+S2*C1 
20  2X  I  =  04  ,  (A ,  04 
J  “=  I  +  J-G4 
01= JO+M 
02=0 l+r 
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R  G  =  >  (03  )  +  X  (  0  2  ) 

F.  1  =  X  (  0  v  )  -  X  (  J  2  ) 
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't  J 


R  3  -  X  (  J  1  )  -  X  (  J  3  ) 
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X(Jv  )  =  P. C +  P? 

Y  (  J  '  )  =  I  C  +  I  2 
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CO  TO  2CC 
ICC  CONTINUE 

X  (  J  2  )  =  P 1  +  I  3 

Y  (  J  2  )  -  I  1  -  R3 
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,0  CONTINUE 
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Y  {  I  )  =  I C 
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6CC  CONTINUE 
7C0  CONTI  NUL 
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DO  30 C  K=2,I 2 
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CC  C  CONTINUE 
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DO  ICC  K - 2  i  13 
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CUM  1 1 .  U  L 
A  L  =  B  S 
J  J  =  0 

GG  2 j C  A  = 1 i A  L 
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DO  20 C  C=8  »  C  L  » C  S 
DC.  20C  i)=C  , DL , D  S 
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2C(  CONTINUE 
RETURN 
END 


